AFRL-VS-HA-TR-2005-1050 


Improved  Surface  Wave  Dispersion 
Models,  Amplitude  Measurements  and 
Azimuth  Estimates 


J.L.  Stevens  M.G.  Eneva 

D.A.  Adams  H.  Xu 

G.E.  Baker 


Science  Applications  International  Corporation 
10260  Campus  Point  Drive 
San  Diego,  CA  92121-1578 


13  March  2005 


Final  Report 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AIR  FORCE  RESEARCH  LABORATORY 
Space  Vehicles  Directorate 
29  Randolph  Rd 

AIR  FORCE  MATERIEL  COMMAND 
Hanscom  AFB,  MA  01731-3010 


This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


AFRL-VS-HA-TR-2005- 1050 


/signed/ 


/signed/ 


ROBERT  A.  RAISTRICK 


ROBERT  BELAND,  Chief 
Battlespace  Surveillance  Innovation  Center 


Contract  Manager 


This  report  has  been  reviewed  by  the  ESC  Public  Affairs  Office  (PA)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS). 

Qualified  requestors  may  obtain  additional  copies  form  the  Defense  Technical 
Information  Center  (DTIC).  All  others  should  apply  to  the  National  Technical 
Information  Service. 

If  your  address  has  changed,  if  you  with  to  be  removed  from  the  mailing  list,  or  if 
the  addressee  is  no  longer  employed  by  your  organization,  please  notify 
AFRL/VSIM,  29  Randolph  Rd„  Hanscom  AFB,  MA  01731-3010.  This  will  assist 
us  in  maintaining  a  current  mailing  list. 

Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on  a 
specific  document  require  that  it  be  returned. 

Using  Government  drawings,  specifications,  or  other  data  included  in  this  document 
for  any  purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the 
U.S.  Government.  The  fact  that  the  Government  formulated  or  supplied  the  drawings, 
specifications,  or  other  data  does  not  license  the  holder  or  any  other  person  or 
corporation;  or  convey  any  rights  or  permission  to  manufacture,  use,  or  sell  any 
patented  invention  that  may  relate  to  them. 

This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange 
and  its  publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its 
ideas  or  findings. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 

data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense.  Washington  Headquarters  Services.  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway.  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  disolav  a  currently 
valid  OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 

31-03-2005 


2.  REPORT  TYPE 

Final 


4.  TITLE  AND  SUBTITLE 

Improved  Surface  Wave  Dispersion  Models,  Amplitude  Measurements  and  Azimuth  Estimates 


6.  AUTHOR(S) 

Jeffry  L.  Stevens,  David  A.  Adams,  G.  Eli  Baker,  Mariana  G.  Eneva,  and  Heming  Xu 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Science  Applications  International  Corporation 
10260  Campus  Point  Drive 
San  Diego,  C A  92121 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited. 


3.  DATES  COVERED  (From  -  To) 

25  Sept  2001  -31  March  2005 


5a.  CONTRACT  NUMBER 

DTRA01-01-C-0082 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 

DTRA 


5e.  TASK  NUMBER 

OT 


5f.  WORK  UNIT  NUMBER 

A1 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/VSBYE 


11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

AFRL-VS-HA-TR-2005- 1 050 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 


The  goal  of  this  project  is  to  reduce  the  magnitude  threshold  for  which  surface  waves  can  be  identified  and 
measured  reliably,  and  to  improve  the  accuracy  of  surface  wave  measurement,  using  phase-matched  filtering  and 
global  regionalized  earth  and  dispersion  models.  Significant  products  and  results  of  this  work  include:  1 )  assembly 
of  a  data  set  of  dispersion  measurements  of  over  one  million  data  points;  2)  tomographic  inversion  of  this  data  set 
for  global  earth  and  dispersion  models;  3)  implementation  and  testing  of  an  improved  azimuth  estimation  technique 
using  Rayleigh  wave  polarization;  4)  implementation  and  testing  of  a  path  corrected  spectral  magnitude.  The  path 
corrected  spectral  magnitude  is  a  “regional  Ms”  which  has  been  a  long-term  goal  of  this  program.  The  global  earth 
and  dispersion  models  are  available  on  request  to  other  researchers  working  in  this  program. 


15.  SUBJECT  TERMS  “  "  - - 

Rayleigh  wave.  Dispersion,  Magnitude,  Polarization,  Global  tomography 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

SAR 

18.  NUMBER 

OF  PAGES 

72 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Robert  J.  Raistrick 

a.  REPORT 

UNCLAS 

b.  ABSTRACT 

UNCLAS 

C.  THIS  PAGE 

UNCLAS 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

781-377-3726 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


ABSTRACT 


The  goal  of  this  project  is  to  reduce  the  magnitude  threshold  for  which  surface  waves  can  be 
identified  and  measured  reliably,  and  to  improve  the  accuracy  of  surface  wave  measurement, 
using  phase-matched  filtering  and  global  regionalized  earth  and  dispersion  models.  Significant 
products  and  results  of  this  work  include:  1)  assembly  of  a  data  set  of  dispersion  measurements 
of  over  one  million  data  points;  2)  tomographic  inversion  of  this  data  set  for  global  earth  and 
dispersion  models;  3)  implementation  and  testing  of  an  improved  azimuth  estimation  technique 
using  Rayleigh  wave  polarization;  4)  implementation  and  testing  of  a  path  corrected  spectral 
magnitude.  The  path  corrected  spectral  magnitude  is  a  “regional  Ms”  which  has  been  a  long-term 
goal  of  this  program.  The  global  earth  and  dispersion  models  are  available  on  request  to  other 
researchers  working  in  this  program. 


TABLE  OF  CONTENTS 

1  Executive  Summary . 1 

2  Overview . 3 

3  Global  Tomographic  Inversion  for  Earth  Models  and  Dispersion  Maps . 5 

3.1  Dispersion  Data  Set . 5 

3.2  Global  Earth  Models  and  Dispersion  Models . 7 

3.3  The  Inversion  Procedure  for  the  3D  Earth  Model . 7 

3.4  Regularization  and  Data  Fit . 8 

3.5  Identification  of  Data  Problems  and  Improvements  to  Data  Quality . 8 

3.6  Prediction  Statistics  for  One-Degree  Earth  Model . 9 

3.7  Surface  Wave  Amplitude  Predictions  from  Global  Earth  Models . 11 

3.8  Properties  of  the  One-degree  Earth  Model . 1 1 

4  Optimized  Surface  Wave  Amplitude  Measurements . 22 

4.1  Path  Corrected  Spectral  Magnitudes . 22 

4.2  Surface  Wave  Measurements  Using  Events  On  and  Near  the  Lop  Nor  Test  Site . 24 

4.3  Path  Corrected  Noise  Estimates  at  the  Lop  Nor  Test  Site . 28 

4.4  Comparison  of  Spectral  and  Time  Domain  Magnitudes . 30 

5  Backazimuth  Estimation  Reliability  Using  Surface  Wave  Polarization . 32 

5.1  Differences  Between  the  Azimuth  Estimation  Algorithms . 32 

5.1.1  The  Current  Method . 32 

5.1.2  The  Chael/Selby  Algorithm  and  its  Implementation . 33 

5.2  Performance  of  the  Algorithms  in  4  Frequency  Bands . 33 

5.3  The  Effect  of  Signal  Strength  (Event  Size)  on  Backazimuth  Estimates . 35 

5.4  The  Effect  of  Variable  Group  Velocity  Windows . 36 

5.5  Predicting  Accuracy  of  Estimates  Using  the  Cross-Correlation  Value . 37 

5.6  Performance  on  Events  for  Which  Surface  Waves  Are  Not  Detected . 38 

6  Performance  of  Algorithms  and  Analysis  of  “Problem  Cases” . 40 

7  Conclusions  and  Recommendations . 51 


v 


8  Data  Deliverable . 52 

9  References . 53 


vi 


LIST  OF  FIGURES 

Figure  1.  Explosions  and  earthquakes  plotted  on  an  Ms:mb  diagram  (from  Stevens  and 
McLaughlin,  2001).  Open  symbols  are  upper  bounds  on  Ms . 3 

Figure  2.  Bar  graph  of  the  number  of  group  velocity  measurements  (left)  and  phase 
velocity  measurements  (right)  in  each  frequency  band  for  data  currently  used  in  the 
tomographic  inversions . 6 

Figure  3.  Distribution  of  group  velocity  data  at  50  seconds  (top  left),  30  seconds  (top 
right),  20  seconds  (lower  left)  and  10  seconds  (lower  right).  The  maps  show  hit  densities 
in  log  10  scale  for  each  one  degree  cell.  Red  indicates  a  high  density  of  coverage  with 
many  rays  crossing  each  cell.  Dark  blue  indicates  little  or  no  data . 6 

Figure  4.  Standard  deviations  (o)  and  means  (+)  of  normalized  group  velocity  residuals 
are  plotted  against  frequency  for  1  degree  earth  model  (solid  red)  and  5-degree  earth 
model  (dashed  blue) . 9 

Figure  5.  Comparison  of  model  misfit  and  data  error.  The  top  figure  shows  a  map  of 
sets  of  closely  spaced  paths,  each  containing  a  bundle  of  at  least  14  rays.  The  b  figure 
shows  data  residuals  relative  to  the  model  (red,  +  are  averages  and  o  standard  deviations) 
and  standard  deviations  of  the  measurements  along  similar  paths  (blue).  Black  squares 
are  the  mean  of  the  standard  deviations  at  the  same  frequency . 10 

Figure  6.  Average  group  velocity  plotted  vs.  frequency  and  ±  one  standard  deviation. 

Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and 
oceanic  models  (blue) . 12 

Figure  7.  Map  of  group  velocities  at  50  seconds  (top)  and  20  seconds  (bottom) . 13 

Figure  8.  Average  phase  velocity  plotted  vs.  frequency  and  ±  one  standard  deviation. 

Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and 
oceanic  models  (blue) . 14 

Figure  9.  Map  of  phase  velocities  at  50  seconds  (top)  and  20  seconds  (bottom) . 15 

Figure  10.  Average  attenuation  coefficient  plotted  vs.  frequency  and  ±  one  standard 
deviation.  Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models 
(red)  and  oceanic  models  (blue) . 16 

Figure  1 1 .  Map  of  phase  velocities  at  50  seconds  (top)  and  20  seconds  (bottom) . 17 

Figure  12.  Average  source  amplitude,  SI,  plotted  vs.  frequency  and  ±  one  standard 
deviation.  Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models 
(red)  and  oceanic  models  (blue) . 1 8 

Figure  13.  Map  of  SI  at  50  seconds  (top)  and  20  seconds  (bottom) . 19 

vii 


Figure  14.  Average  source  amplitude,  S2,  plotted  vs.  frequency  and  ±  one  standard 
deviation.  Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models 
(red)  and  oceanic  models  (blue) . 20 

Figure  1 5.  Map  of  S2  at  50  seconds  (top)  and  20  seconds  (bottom) . 21 

Figure  16.  Path  corrected  spectra  for  an  explosion  and  for  earthquakes  calculated  for 
several  depths.  The  path  corrected  explosion  spectrum  is  flat  over  the  entire  frequency 
band  (for  perfect  data  and  path  correction),  while  the  path  corrected  earthquake  spectrum 
is  flattened,  but  has  some  variation  due  to  source  mechanism  and  source  depth . 23 

Figure  17.  Path  corrected  spectral  magnitude  (Log  M0)  residuals  plotted  vs.  distance 
(from  Stevens  and  McLaughlin,  2001).  Log  Mo  is  nearly  independent  of  distance . 23 

Figure  18.  Maps  showing  the  Lop  Nor  area  (rectangle),  stations  (triangles),  and 
earthquake  (circles)  and  explosion  (crosses)  epicenters . 24 

Figure  19.  Examples  of  path  corrected  spectra  used  in  this  work:  (a)  Lop  Nor  explosions 
recorded  at  distances  of  2°,  7°  and  67°  (left);  (b)  Lop  Nor  earthquakes  recorded  at  distances 
of  0.4°,  22°  and  65°  (right).  See  examples  of  station  logMt  estimates  in  Table  1 .  S/N  is  good 
at  all  but  the  highest  frequencies.  Low  frequency  noise  dominates  over  the  surface  wave 
signal  below  about  0.02  Hz . 25 

Figure  20.  Comparison  of  station  spectral  magnitudes  calculated  with  two  different 
methods.  Insets  show  standard  deviations  as  indicated.  See  text  for  details . 26 

Figure  21.  Histograms  of  standard  deviations  of  the  mean  station  logAT0  estimated  in 
three  frequency  bands  (shown  on  right),  for  small  and  larger  distances  (shown  on  top) . 26 

Figure  22.  Comparison  of  station  spectral  magnitudes  in  different  frequency  bands,  (a) 
adjacent  bands,  all  distances;  (b)  overlapping  bands,  all  distances;  (c)  overlapping  bands, 
distances  <5° . . . 27 

Figure  23.  Log/V/o:/Mb  plots  showing  event  spectral  magnitudes  for  earthquakes  (O)  and 
explosions  (X)  in  Lop  Nor.  Station  spectral  magnitudes  were  calculated  using  frequencies 
increasing  with  decreasing  distance  (left)  and  the  0.03-0.07  Hz  frequency  band  for  all 
distances  (right).  Bold  lines  mark  the  empirical  discrimination  relationship  of  Stevens  and 
McLaughlin  (2001).  For  small  events  higher  frequencies  may  provide  the  best  signal  to 
noise  performance,  but  at  a  cost  of  degraded  discrimination  capability . 28 

Figure  24.  Average  “path  corrected  noise  measurement”  and  ±one  standard  deviation 
curves  for  13  time  segments  at  WMQ  (left)  and  for  54  time  segments  at  HI  A  (right).  The 
path  corrected  noise  measurement  gives  the  minimum  path  corrected  signal  level  that 
could  be  observed  at  each  frequency.  Recall  that  the  predicted  path  corrected  explosion 
spectrum  is  flat  over  this  frequency  band,  and  the  predicted  path  corrected  earthquake 
spectra  decrease  with  increasing  frequency  (see  Figure  16) . 29 


viii 


Figure  25.  Backazimuth  residuals  for  0.03-0.05  Hz  using  the  current  method  (left)  and 
the  Chael/  Selby  algorithm  using  the  normalization  of  equation  5.2  and  peak  cross- 


correlation  (right) . 34 

Figure  26.  Median  ±2  SMAD  confidence  intervals  of  backazimuth  residuals  binned  by 
M,  values  for  the  CM  (dotted)  and  CSi  algorithm  (solid).  Each  bin  has  250  backazimuth 
residuals.  Results  at  0.03-0.05  Hz  are  intermediate  between  those  shown . 36 


Figure  27.  Surface  waves  in  the  3  passbands  tested,  at  the  station  ARCES,  for  an  Ms  5.0 
event  at  6035  km.  The  entire  trace  is  the  5.0  to  2.5  km/s  group  velocity  window.  Outlined 
segments  are  the  shorter  group  velocity  windows  used  to  isolate  the  Rayleigh  waves.  The 
isolation  works  best  in  the  highest  frequency  band,  where  the  large  Love  wave  (top  trace)  is 


outside  the  narrower  window . 37 

Figure  28  Median  azimuth  residuals  ±2  SMAD  confidence  intervals  vs.  the  correlation 
of  the  radial  and  Hilbert  transformed  vertical  Rayleigh  waves  for  the  0.03-0.05  Hz 
passband.  Results  are  similar  for  the  other  passbands . 38 


Figure  29.  Azimuth  residuals  for  the  0.02  -  0.04  Hz  passband  for  the  CS  algorithm  (top 
row)  and  the  CM  (bottom  row).  Histograms  of  azimuth  residuals  for  all  the  data  (left 
column)  and  for  data  with  correlation  £  0.8  for  CS  (upper  middle)  or  F-statistic  ^  30  for 
the  CM  (lower  middle).  The  middle  plots  use  -14.5%  of  the  data  in  each  case.  The  upper 
right  plot  shows  the  median  azimuth  residual  ±2  SMAD  confidence  intervals  vs.  the 


correlation,  as  in  Figure  28,  for  the  CS  algorithm.  The  bottom  right  shows  the  median 
azimuth  residual  ±2  SMAD  confidence  vs.  the  F  statistic  of  the  CM . 39 

Figure  30.  Seismograms  recorded  at  two  stations  on  paths  passing  through  the  Tarim 
Basin.  There  are  two  distinct  surface  wave  arrivals  of  comparable  amplitude.  The  second 
arrival  is  dominant  for  paths  traveling  directly  through  the  basin,  and  the  first  arrival  is 
dominant  for  paths  near  the  basin  boundary . 40 

Figure  31.  Paths  near  and  across  the  Tarim  Basin.  Colors  correspond  to  the  colors  on  the 
dispersion  curves  in  Figure  32 . 4 1 

Figure  32.  Observed  (left)  and  calculated  (right)  dispersion  curves  for  the  paths  shown  in 
Figure  3 1 .  Colors  correspond  to  the  paths  on  the  map  in  Figure  31 . 41 

Figure  33.  3-component  surface  waves  for  recordings  of  a  single  event  with  one  path 
crossing  (right)  and  one  not  crossing  (left),  the  Tarim  Basin.  Source-station  distances  are 
1,427  km  (left)  and  1,918  km  (right) . 42 

Figure  34.  Propagation  paths  overlain  on  a  sediment  thickness  map  (contoured  in 
kilometers  by  grayscale  darkness)  for  the  uncomplicated  records  on  the  left  side  of  Figure 
6.4  (propagating  to  the  northeast)  and  the  complicated  records,  which  propagate  west 
through  the  thick  sediments  of  the  Tarim  Basin . 43 

Figure  35.  Dispersion  curves  (red)  picked  for  the  vertical  component  records  shown  in 
Figure  33.,  compared  with  the  model  based  prediction  (maroon) . 44 


Figure  36.  Vertical  component  surface  waves  (top  row)  for  the  problem  path  across  the 
Tarim  Basin.  The  second  row  shows  the  correlation  of  the  Hilbert  transformed  radial  with 
the  vertical,  which  identifies  Rayleigh  waves.  The  third  row  shows  the  backazimuth  vs 
time.  The  left  column  is  for  data  filtered  from  10  to  20  seconds  period,  and  the  right  is  for 


data  filtered  from  5  to  10  seconds  period.  The  great  circle  back  azimuth  is  101  degrees . 45 

Figure  37.  Location  of  paths  through  different  parts  of  the  Tarim  Basin  (red  lines)  from 
events  (within  green  circle  at  lower  right)  to  the  KNET  station  ULHL  (blue  circle). 
Sediment  thickness  in  kilometers  is  indicated  by  grayscale  shading . 46 


Figure  38.  3-component  records,  filtered  from  20  to  30  second  period,  arranged  by 
backazimuth  of  the  events  within  the  green  circle  in  the  lower  right  of  Figure  6.8 
Records  are  centered  at  2.8  km/sec  (0  seconds  on  the  time  axis).  Fifty  seconds  either  side 


of  the  center  are  highlighted  in  each  record  to  facilitate  comparison . 47 

Figure  39.  Same  as  Figure  38,  except  filtered  from  10  to  20  second  period . 47 

Figure  40.  Vertical  (blue)  and  Hilbert  transformed  radial  (dotted  red)  seismograms 
(left)  filtered  from  20  to  30  seconds  period,  and  their  correlation  over  a  series  of  tightly 
overlapping  75  second  windows . 48 

Figure  41.  Vertical  (blue)  and  Hilbert  transformed  radial  (dotted  red)  seismograms 
(left)  filtered  from  10  to  20  seconds  period,  and  their  correlation  over  a  series  of  tightly 
overlapping  45  second  windows . 49 


Figure  42.  Vertical  component  surface  waves  (top  row)  for  the  two  southernmost 
records  of  Figure  6.12  (126.9°  backazimuth  on  the  left,  125.0°  backazimuth  on  the  right), 
the  correlation  of  the  horizontal  component  (Hilbert  transformed)  rotated  to  the 
backazimuth  that  maximizes  the  crosscorrelation  with  the  vertical  for  successive 
windows  (middle  row),  and  that  backazimuth  vs  time  (bottom  row).  The  horizontal  line 
indicates  the  theoretical  backazimuth.  Data  are  filtered  from  10  to  20  seconds  period. 

Time  is  seconds  from  origin . 50 


x 


LIST  OF  TABLES 


Page 

Table  1 .  Station  LogMO  estimates  in  0.05-0. 10  Hz  from  the  spectra  in  Figure  19 . 27 

Table  2.  Comparison  of  time  domain  and  spectral  magnitude  measurement  and  correction 
terms . 30 

Table  3.  Scaled  median  absolute  deviations  and  standard  deviations  (in  parentheses)  of 
the  backazimuth  residuals  in  4  frequency  bands  from  the  algorithms  tested.  Four 
variations  of  the  implementation  of  the  CS  algorithm  are  tested.  CSi  and  CS2  use  the  peak 
amplitude  of  the  cross-correlation  while  CS?  and  CS4  use  the  circular  mean.  CSi  and  CS? 

are  normalized  by  -Js^S~  as  in  equation  5.2.  CS2  and  CS4  are  normalized  by  S!Z  as  in 
equation  5.4 . 35 

Table  4.  Same  as  Table  3,  but  for  group  velocity  windows  designed  to  isolate  the 
Rayleigh  waves.  The  change  from  the  2.5  to  5  km/s  group  velocity  windows  is  given  in 
parentheses  (negative  change  is  improvement) . 37 


xi 


PREFACE 


We  would  like  to  thank  Mike  Ritzwoller,  Anatoli  Levshin  and  Nikolai  Shapiro  of  the  University 
of  Colorado,  Bob  Herrmann  of  St.  Louis  University,  and  other  researchers  who  have  allowed  us 
to  use  their  data  in  this  project.  Much  of  the  work  described  in  this  report  has  been  presented  in 
three  Seismic  Research  Review  papers  (Stevens  et  al,  2002,  2003  and  2004).  The  material  in 
section  6  has  been  published  in  Geophysical  Research  Letters  (Baker  and  Stevens,  2004). 


1  EXECUTIVE  SUMMARY 

The  primary  goal  of  this  project  is  to  reduce  the  magnitude  threshold  for  which  surface  waves 
can  be  identified  and  measured  reliably,  and  to  improve  the  accuracy  of  surface  wave 
measurement,  using  phase-matched  filtering  and  global  regionalized  earth  and  dispersion 
models.  Following  are  the  most  significant  products  and  results  of  this  work: 

1 .  We  assembled  a  dataset  of  dispersion  measurements  of  over  one  million  data  points. 
Some  of  this  dataset  came  from  our  own  measurements,  but  most  came  from  other 
research  studies,  and  we  thank  all  of  those  scientists  who  have  contributed  to  this  effort. 
This  dataset  is  one  of  the  largest  dispersion  datasets  ever  assembled  and  is  particularly 
unique  in  its  coverage  for  a  broad  range  of  frequencies. 

2.  We  developed  and  have  continuously  improved  a  set  of  global  earth  models  and 
dispersion  models  defined  on  a  one-degree  grid.  These  models  were  developed  by 
simultaneously  inverting  the  entire  dataset  for  a  set  of  earth  structures,  which  in  turn 
allow  surface  wave  dispersion  to  be  calculated  between  any  two  points  on  the  earth  at  any 
set  of  frequencies,  with  the  dispersion  models  constrained  by  the  large  dataset.  These 
global  earth  and  dispersion  models  are  available  on  request  to  other  researchers  working 
in  this  program. 

3.  We  implemented  the  azimuth  estimation  technique  developed  by  Selby  using  Chad's 
algorithm,  tested  it  on  a  large  dataset  of  International  Monitoring  System  (IMS)  data,  and 
demonstrated  that  it  is  a  major  improvement  over  other  commonly  used  azimuth 
estimation  techniques.  This  azimuth  estimation  technique  can  be  used  together  with  the 
existing  detection  test,  which  is  based  on  consistency  of  measurement  with  surface  wave 
dispersion  maps,  to  reduce  the  detection  threshold  for  which  surface  waves  can  be 
reliably  identified  and  measured. 

4.  We  implemented  and  tested  a  path  corrected  spectral  magnitude,  developed  procedures 
for  optimizing  it  and  evaluated  the  discrimination  capability  of  this  and  other  types  of 
surface  wave  measurements.  Two  important  results  of  this  work  are: 

a.  The  path  corrected  spectral  magnitude  is  independent  of  distance  and  unaffected  by 
dispersion  and  therefore  can  be  measured  at  any  distance,  including  very  close  to  the 
source,  and  it  will  have  the  same  value  as  a  measurement  made  at  a  greater  distance. 
The  path  corrected  spectral  magnitude  is  therefore  a  “regional  M«”  which  has  been  a 
long-term  goal  of  this  program. 

b.  For  discrimination  purposes,  surface  waves  should  be  measured  at  periods  greater 
than  1 0  seconds.  Periods  of  1 0  seconds  and  longer  can  be  measured  even  at  very 
close  range,  and  there  is  little  or  no  S/N  improvement  at  shorter  periods. 
Discrimination  is  degraded  at  higher  frequencies  due  to  the  decrease  in  earthquake 
spectral  amplitude  relative  to  explosions  at  higher  frequencies. 
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5.  We  performed  a  study  of  “problem  cases”,  where  surface  wave  dispersion  and/or 
amplitudes  are  inconsistent  with  model  predictions.  Typically  these  cases  occur  where 
there  are  strong  heterogeneities  in  earth  structure  along  the  path,  particularly  for  grazing 
incidence  along  large  changes  in  material  properties.  Ray  theory  predicts  that  surface 
waves  will  take  the  minimum  time  path  between  the  source  and  receiver  in  such  cases, 
but  the  observed  waveforms  are  considerably  more  complex  than  this.  In  cases  of  large 
velocity  contrast  such  as  the  Tarim  Basin  region,  there  are  multiple  surface  wave  arrivals 
with  a  distinct  arrival  passing  through  the  basin  and  another  traveling  around  it.  The 
arrivals  merge  at  longer  periods.  Special  care  is  required  for  construction  and  use  of 
dispersion  and  attenuation  models  in  such  cases. 
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2  OVERVIEW 


The  importance  of  long  period  (>~10  second)  surface  waves  for  nuclear  monitoring  is  in  the 
discrimination  capability  of  the  Ms:mb  discriminant  and  its  variants  (Marshall  and  Basham,  1972; 
Stevens  and  Day,  1985).  In  general,  contained  underground  explosions  generate  larger  mb 
relative  to  Ms  than  earthquakes  of  comparable  size  (Figure  1 ). 

Ms:mb 


Figure  I.  Explosions  and  earthquakes  plotted  on  an  Ms:mb  diagram  (from  Stevens  and 
McLaughlin,  2001).  Open  symbols  are  upper  bounds  on  Ms. 

Ms:mb  is  one  of  the  most  reliable  earthquake/explosion  discriminants,  but  it  has  some  limitations 
which  we  have  tried  to  address  in  this  study: 

1 .  Body  waves  can  be  measured  for  significantly  smaller  events  than  surface  waves,  so  the 
discriminant  is  not  useful  for  very  small  events  We  have  addressed  this  issue  by 
investigating  methods  to  reduce  the  threshold  for  which  surface  waves  can  be  reliably 
identified  and  measured,  including  development  of  improved  dispersion  maps 
constrained  by  a  very  large  data  set  to  facilitate  phase-matched  filtering. 

2  While  the  strongest  surface  wave  signals  are  at  the  closest  distances,  Ms  is  an  unreliable 
amplitude  measure  at  regional  and  local  distances  because  of  differences  in  dispersion 
and  variations  in  the  frequency  content  of  the  signal.  We  have  addressed  this  issue  by 
developing  surface  wave  measurement  techniques  that  can  be  used  at  local  and  regional 
distances,  even  very  close  to  the  source,  which  are  consistent  with  measurements  made  at 
greater  distances. 
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3.  The  discrimination  capability  of  surface  waves  decreases  with  increasing  frequency 
because  the  spectral  difference  between  earthquakes  and  explosions  declines.  We  have 
addressed  this  issue  by  analyzing  surface  wave  signal  and  noise  measurements  to  define 
the  optimum  frequency  bands  for  discrimination. 

4.  Association  of  surface  waves  with  the  wrong  event  is  a  problem  that  can  potentially  lead 
to  misclassification  of  an  explosion  as  an  earthquake.  We  have  addressed  this  issue  in 
two  ways:  first,  by  improving  dispersion  models  that  are  used  in  the  dispersion  test  for 
surface  wave  identification;  and  second  by  implementing  and  testing  an  improved 
azimuth  estimation  technique  that  can  be  used  to  improve  association  of  surface  waves 
with  the  corresponding  event. 

5.  In  general  surface  waves  can  be  modeled  quite  well  using  path-averaged  dispersion  and 
attenuation  calculated  from  discrete  plane-layered  earth  models  and  great  circle 
propagation,  however  there  are  some  complex  cases  where  standard  methods  of 
measuring  surface  wave  dispersion  and  amplitude  measurement  do  not  work  well  We 
have  investigated  some  problem  cases  and  made  recommendations  for  techniques  to 
handle  these  situations. 
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3  GLOBAL  TOMOGRAPHIC  INVERSION  FOR  EARTH  MODELS 
AND  DISPERSION  MAPS 

We  develop  global  regionalized  earth  and  dispersion  models  by  inversion  of  a  very  large  data  set 
of  phase  and  group  velocity  dispersion  measurements.  The  data  set  has  grown  from  about  90,000 
(Stevens  and  McLaughlin,  2001)  to  over  one  million  data  points,  and  the  resolution  has  improved 
such  that  the  base  model  has  changed  from  a  5  degree  grid  to  a  one  degree  grid.  The  dispersion 
measurements  are  inverted  for  earth  structure,  and  the  earth  structure  is  then  used  to  generate 
dispersion  predictions  as  described  below.  The  phase  velocity  dispersion  curves  are  then  used  to 
calculate  phase-matched  filters  to  improve  detection.  An  important  advantage  gained  from  using 
earth  models  is  that  we  can  include  information  from  other  studies  leading  to  physically 
reasonable  constraints  on  dispersion.  For  our  earth  models  this  information  consists  of  the 
boundaries  between  geologic  zones,  bathymetry  of  oceans,  thicknesses  of  sediments  and  ice, 
Moho  depths,  and  prior  estimates  of  seismic  velocities  derived  from  Crust  2.0  and  AK  135  earth 
models.  These  constraints  are  important  for  filling  in  the  gaps  found  in  the  path  coverage  of  our 
data  set  and  they  enable  dispersion  prediction  along  paths  with  limited  data. 

3.1  Dispersion  Data  Set 

The  dispersion  data  set  has  been  derived  from  a  variety  of  regional  and  global  studies  including  the 
following:  global  surface  wave  group  velocities  from  earthquakes  derived  using  IMS  data  (Stevens 
and  McLaughlin,  1996),  augmented  with  additional  measurements  derived  from  IMS  data;  surface 
wave  phase  and  group  velocity  dispersion  curves  from  underground  nuclear  test  sites  (Stevens, 
1986;  Stevens  and  McLaughlin,  1988),  calculated  from  earth  models  for  270  paths  (test  site  - 
station  combinations)  at  10  frequencies  between  0.015  and  0.06  Hz;  phase  and  group  velocity 
measurements  for  western  Asia  and  Saudi  Arabia  from  Mitchell  et  al.(  1 996)  for  12  paths  at  17 
frequencies  between  0.012  and  0. 14  Hz;  the  global  phase  velocity  model  of  Ekstrom  et  al.  (1996) 
for  9  periods  between  35  and  1 50  seconds  calculated  for  each  grid  block  from  a  spherical  harmonic 
expansion  of  order  1  =  40;  group  velocity  measurements  for  Eurasia  from  Ritzwoller  et  al.  (1996) 
and  Levshin  et  al  (1996)  for  20  frequencies  between  0.004  and  0. 1  Hz  with  500  to  5000  paths  per 
frequency;  Antarctic  and  South  American  group  velocity  measurements  from  the  University  of 
Colorado  (Vdovin  et  al.,  1999;  Ritzwoller  et  al.,  1999);  high  frequency  Eurasian  dispersion 
measurements  from  University  of  Colorado  (Levshin  et  al,  2003);  dispersion  measurements  from 
Central  Asia  made  by  Los  Alamos  National  Laboratory  (Yang  et  al,  2002);  data  from  China 
(Huang  et  al,  2003);  and  data  from  the  New  Madrid  region  from  Mancilla  (2001). 

Figure  2  shows  the  frequency  distribution  of  group  velocity  and  phase  measurements  in  our 
data  set,  excluding  the  phase  velocities  derived  from  the  global  phase  velocity  model  of  Ekstrom 
et  al.  (1996).  Figure  3  shows  the  geographic  distribution  of  data  coverage.  Coverage  is 
excellent  in  Eurasia  for  all  frequencies,  and  good  in  most  continental  regions  at  20  seconds  and 
lower.  Data  coverage  is  more  limited  in  the  oceans,  particularly  at  higher  frequencies.  The  maps 
show  “hit  densities,”  in  loglO  scale,  over  the  one  degree  cells  for  all  group  velocity 
measurements  in  the  data  for  the  indicated  frequency  ranges.  The  density  for  each  cell  is  the  sum 
of  proportions  of  the  lengths  to  total  lengths  of  all  paths  passing  through  the  cell.  So,  for 
example,  a  hit  density  of  10  could  mean  that  100  rays,  each  10  degrees  in  length,  pass  through 
the  cell. 
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Figure  2.  Bar  graph  of  the  number  of  group  velocity  measurements  (left)  and  phase  velocity 
measurements  (right)  in  each  frequency  band  for  data  currently  used  in  the  tomographic 
inversions. 


0.0453-0.05 17  Hz  0.1054 -0.1 15  Hz 

Figure  3.  Distribution  of  group  velocity  data  at  50  seconds  (top  left),  30  seconds  (top  right),  20 
seconds  (lower  left)  and  10  seconds  (lower  right).  The  maps  show  hit  densities  in  log  10  scale  for 
each  one  degree  cell.  Red  indicates  a  high  density  of  coverage  with  many  rays  crossing  each  cell. 
Dark  blue  indicates  little  or  no  data. 
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3.2  Global  Earth  Models  and  Dispersion  Models 

The  regionalized  earth  model  consists  of  l°xl°  blocks  and  is  made  up  of  layers  of  ice,  water, 
sediments,  crust  and  upper  mantle.  This  model  depends  on  9022  free  parameters  which  are 
adjusted  by  a  damped  least  squares  fit  to  Rayleigh  wave  dispersion  data.  The  free  parameters  are 
the  S-wave  velocities  of  layers  of  577  different  model  types.  Other  constrained  parameters  in  the 
model  are  P  wave  velocities,  densities  and  Q.  The  model  types  are  based  on  the  Crust  2.0  2°x2° 
crustal  types  (Bassin  et  al.,  2000  and  Laske  et  al.  2001)  and  also  on  ocean  ages  (Stevens  and 
Adams,  2000).  The  top  few  km  of  the  model  (consisting  of  water,  ice  and/or  sediments)  are  fixed 
and  match  data  from  one  degree  bathymetry  maps  made  by  averaging  Etopo5  five  minute 
measurements  of  topography,  and  Laske  and  Masters  ( 1 997)  1  degree  maps  of  sediments.  There 
is  an  explicit  discontinuity  between  the  bottom  of  the  sediments  and  the  crust.  There  are  three  or 
more  crustal  layers.  The  Crust  2.0  models,  which  were  the  starting  point  for  these  structures, 
have  three  crustal  layers,  but  we  found  it  necessary  to  add  more  layers  in  regions  of  thick  crust. 
There  is  another  explicit  discontinuity  at  the  Crust/Mantle  boundary.  The  Moho  depth  is  derived 
from  Crust  2.0  and  varies  on  a  2°  grid.  The  mantle  starting  model  is  derived  from  AK135 
(Kennett,  et  al,  1995).  With  these  constraints,  the  inversion  is  performed  for  the  shear  velocity  of 
the  crust  and  upper  mantle  to  a  depth  of  3 1 0  km.  Below  310km  the  earth  structure  is  fixed,  and 
the  inversion  model  is  required  to  be  continuous  with  the  mantle  structure  at  the  base  of  the 
inversion.  In  broad  ocean  areas,  we  replace  the  Crust  2.0  model  with  models  distinct  for  each 
ocean  and  subdivided  by  ocean  age.  We  also  separate  into  distinct  models  Crust  2.0  models  that 
are  geographically  separated.  So,  for  example,  if  Crust  2.0  has  the  same  model  type  in  North 
America  and  in  Asia,  we  use  the  same  starting  model  for  each,  but  treat  them  as  separate  models 
in  the  inversion. 

3.3  The  Inversion  Procedure  for  the  3D  Earth  Model 

The  relationship  between  dispersion  and  the  shear  wave  velocities  of  the  layers  in  the  earth 
model  is  non-linear,  so  the  shear  velocities  are  estimated  by  an  iterative  least  squares  inversion 
procedure.  At  each  step  a  system  of  tomographic  equations  is  formed,  augmented  by  additional 
equations  of  constraint,  and  then  solved  by  the  LSQR  algorithm.  The  equations  solved  are 


fA' 

f  Ad  ' 

sH 

Ax  = 

-sHx 

UJ 

where  Ax  is  the  vector  of  adjustments  to  the  shear  wave  slownesses  of  layers  in  each  of  the  577 
model  types.  Ad  is  the  vector  of  slowness  differences  between  predicted  and  observed  dispersion 
measurements,  e  is  the  vector  of  residuals  that  remain  after  inversion  (the  inversion  minimizes 
|e|-).  x  is  the  vector  of  slownesses  estimated  in  the  previous  iteration.  The  elements  of  the  matrix 
A  consist  of  partial  derivatives  of  dispersion  predictions  with  respect  to  shear  wave  slownesses  in 
each  layer.  H  is  a  difference  operator  that  applies  to  vertically  neighboring  layers  and  has  the 
effect  of  constraining  the  vertical  smoothness  of  velocity  profile.  H  applies  to  layers  in  the  crust 
and  upper  mantle,  but  has  explicit  discontinuities  at  the  crust/mantle  boundary  and  at  the  base  of 
surface  sediments,  s  is  the  weighting  of  the  smoothness  constraint  and  can  be  a  diagonal  matrix 
(for  variably  weighted  smoothing)  or  a  scalar  (constant  smoothing).  A  different  smoothing 
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parameter  can  be  selected  for  each  model  type.  Lateral  smoothing,  which  is  usually  applied  in 
tomography  studies,  is  executed  indirectly  in  our  study  through  selection  of  the  model  types.  I  is 
the  identity  matrix  and  X  weights  the  damping  which  constrains  the  norm  of  the  difference 
between  final  slownesses  and  constraining  model  slownesses  xc  (in  this  case  a  variant  of  the 
Crust  2.0  values).  X  can  be  a  scalar  for  constant  damping,  or  a  diagonal  matrix  for  variable 
damping.  As  with  smoothing,  variable  damping  is  implemented  so  that  a  different  parameter  can 
be  selected  for  each  model  type. 

An  important  feature  of  Equation  3.1,  and  what  makes  this  procedure  different  from  most  other 
tomographic  inversion  studies,  is  that  we  invert  all  relevant  data  in  the  same  inversion.  The  data 
includes  phase  velocity  and  group  velocity  measurements  measured  along  specific  paths  at  all 
frequencies,  as  well  as  some  phase  and  group  velocity  data  points  derived  from  models  (e  g.  the 
Harvard  phase  velocity  models).  It  is  also  possible  to  include  earth  models  as  data,  so  that  earth 
models  from  existing  studies  can  be  included  in  the  inversion. 

3.4  Regularization  and  Data  Fit 

The  damping  and  smoothing  constraints  and  their  associated  weighting  parameters  are  used  to 
regularize  the  solution.  Regularization  acts  both  to  control  the  influence  of  data  noise  on  the 
estimation  of  model  parameters  and  to  constrain  parts  of  the  model  that  are  poorly  constrained  by 
data.  Too  much  regularization  will  make  the  model  too  smooth,  degrading  the  data  fit,  and  too 
little  regularization  will  allow  noise  to  be  projected  into  the  model,  making  it  unrealistic  and 
rough.  Techniques  for  optimization  of  regularization  parameters  are  not  yet  mature,  especially 
for  large-scale  problems  such  as  this.  The  methods  most  often  used  (e.g.  Hansen,  1998)  are  the 
L-curve,  generalized  cross  validation,  and  discrepancy  principle.  We  have  experimented  with 
several  of  these  techniques  for  our  inversion  problem,  but  have  not  found  any  reliable  enough  to 
replace  analyst  review  of  the  inversion  results. 

3.5  Identification  of  Data  Problems  and  Improvements  to  Data  Quality 

The  quality  of  the  inversion  results  is  limited  by  the  quality  of  the  data.  Redundancy  in  the  data 
helps  to  average  out  errors;  but  as  the  models  improve,  the  remaining  errors  in  the  data  begin  to 
limit  the  improvement  that  is  possible.  Consequently,  we  initiated  a  review  and  assessment  of  the 
data  quality,  and  made  an  effort  to  identify  and  remove  data  problems.  This  has  made  dramatic 
differences  in  the  results  in  some  areas.  Some  of  the  types  of  problems  that  we  identified  are: 

1  Consistent  errors  made  by  research  groups.  In  reviewing  data  sets,  we  find  that  they  often 
contain  consistent  errors  throughout  the  data  set,  most  commonly  errors  in  dispersion  at 
either  the  high  frequency  or  low  frequency  end  of  the  dispersion  curves.  These  errors  can 
be  identified  by  comparison  with  dispersion  curves  from  other  sources,  and  the  problem 
data  range  can  be  removed. 

2.  Errors  in  location.  Particularly  on  short  paths,  event  location  errors  can  lead  to  large 
errors  in  dispersion  curves.  If  we  have  information  about  the  event  used  to  derive  the 
dispersion  curve,  and  a  new,  better  location,  then  we  can  correct  the  error.  We  found,  for 
example,  that  recalculating  group  velocities  using  improved  locations  from  Engdahl  et  al 
(1998)  greatly  improved  the  consistency  of  the  data. 
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3.  Grazing  paths  along  low  velocity  zones.  This  is  a  particularly  troubling  problem,  because 
surface  waves  along  these  paths  often  show  two  distinct  or  interfering  arrivals,  one  that 
travels  through  a  low  velocity  basin  and  one  that  travels  around  it.  Analysts  may  pick  one 
or  the  other  inconsistently.  The  problem  can  be  reduced  by  ordering  the  arrivals  such  that 
the  dispersion  picks  change  from  fast  to  slow  to  fast  smoothly  as  the  azimuth  changes 
across  the  basin  (see  Section  5  of  this  report). 

3.6  Prediction  Statistics  for  One-Degree  Earth  Model 

We  evaluate  the  performance  of  the  one-degree  earth  model  for  predicting  dispersion  in  two 
ways:  1)  by  calculating  the  mean  and  standard  deviation  of  the  data  residual  for  the  entire  data 
set,  and  2)  by  comparing  the  residual  error  in  the  prediction  vs.  data  with  the  consistency  of  the 
data  set  itself. 

The  means  and  standard  deviations  of  normalized  group  velocity  residuals,  1 -Vo/vp,  where  vcand 
vp  are  observed  and  model  predicted  group  velocities,  were  calculated  for  narrow  frequency 
bands  and  are  shown  in  Figure  3.3.  The  solid  line  is  for  our  best  1 -degree  model,  and  dashed  is 
for  the  best  5-degree  model  (Stevens  and  Adams,  2000)  based  on  Crust  5. 1  (Mooney,  et  al., 

1998).  Figure  4  shows  the  value  of  the  1 -degree  model,  especially  for  high  frequencies.  We 
note,  however,  that  the  five  degree  model  was  based  on  a  smaller  data  set  and  the  performance 
would  be  better  if  a  new  inversion  was  performed  with  the  complete,  current  data  set. 


Frequency  (Hz) 

Figure  4.  Standard  deviations  (o)  and  means  (+)  of  normalized  group  velocity  residuals  are 
plotted  against  frequency  for  1  degree  earth  model  (solid  red)  and  5-degree  earth  model  (dashed 
blue). 
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As  earth  models  have  improved  through  the  course  of  this  project,  improvement  in  models  has 
become  more  limited  by  the  accuracy  of  the  data.  Figure  5  illustrates  estimates  of  data  error 
together  with  residuals  from  the  model.  Sets  of  “master”  paths  were  selected  that  each  had  14  or 
more  rays  from  different  events  and/or  stations  along  very  similar  paths.  Each  path  contains  a 
bundle  of  rays  within  a  “distance”  of  .01  of  a  “master”  path.  The  distance  is  (dl+d2)/D  where  D 
is  the  length  of  the  master  path,  and  dl  and  d2  are  the  distances  between  pairs  of  end  points. 
These  figures  show  that  the  data  residuals  are  comparable  to  the  errors  in  the  data,  and  that 
further  improvement  in  data  quality  has  the  potential  to  reduce  data  errors  and  allow  further 
improvements  in  model  prediction. 


Figure  5.  Comparison  of  model  misfit  and  data  error.  The  top  figure  shows  a  map  of  sets  of 
closely  spaced  paths,  each  containing  a  bundle  of  at  least  14  rays.  The  b  figure  shows  data 
residuals  relative  to  the  model  (red,  +  are  averages  and  o  standard  deviations)  and  standard 
deviations  of  the  measurements  along  similar  paths  (blue).  Black  squares  are  the  mean  of  the 
standard  deviations  at  the  same  frequency. 
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3.7  Surface  Wave  Amplitude  Predictions  from  Global  Earth  Models 


Surface  wave  amplitudes  can  be  predicted  with  an  approximation  originally  due  to  McGarr 
(1969)  that  uses  propagation  of  surface  waves  along  great  circle  paths  with  conservation  of 
energy  across  material  interfaces  and  no  mode  conversion.  With  these  approximations,  surface 
wave  propagation  in  a  heterogeneous,  anelastic  structure  takes  the  following  form,  separating 
source,  path  and  receiver  (notation  follows  Harkrider  et  al,  1994): 


u2(<o,r, <p )  =  .  .  1  exp [i ( jt/4 -ar/c  -//)] F, (», <P, h)  (3-2) 

yja,  sm{r/a,)  \  ncoc\  y  L  J 


where  a>  is  angular  frequency,  r  is  source  to  receiver  distance,  ae  is  the  radius  if  the  earth,  ^>is 
azimuth,  Ar  is  the  Rayleigh  wave  excitation  function,  c  is  phase  velocity,  y  is  the  attenuation 
coefficient,  and  the  subscripts  1,  2,  and  p  refer  to  parameters  derived  from  the  source  region 
structure,  parameters  derived  from  the  receiver  region  structure,  and  parameters  which  are 
defined  by  path  averages,  respectively.  All  source  properties  are  contained  in  the  function  Fs.  For 
an  isotropic  explosion  source,  the  Rayleigh  wave  spectrum  can  be  written: 


u,(a>,hx,r)  =  M0 


S?((D,hx)S2(ci))exp[-yp(6))r  +  i(<pt)-cor/cp(o))\ 
yjae  sin(r/tf,) 


(3-3) 


where  (pn  is  the  initial  phase  equal  to  -3rt/4,  S*  depends  on  the  source  region  elastic  structure  and 
the  explosion  source  depth,  hx  depends  on  the  receiver  region  elastic  structure.  M0  =  yA/0 

where  M0  is  the  explosion  isotropic  moment.  Mu  is  defined  this  way  so  that  the  function  S* 

does  not  depend  explicitly  on  the  material  properties  at  the  source  depth  (More  details  are  given 
in  Stevens  and  McLaughlin  (2001)  and  Stevens  and  Murphy  (2001)). 

3.8  Properties  of  the  One-degree  Earth  Model 

In  this  section  we  show  the  values  and  distribution  of  quantities  derived  from  the  one-degree 
earth  model.  Figure  6  shows  group  velocity  as  a  function  of  frequency  for  all  models,  and 
Figure  7  shows  a  map  of  group  velocity  at  50  and  20  seconds. 
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Figure  6.  Average  group  velocity  plotted  vs.  frequency  and  ±  one  standard  deviation.  Top 
figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and  oceanic  models 
(blue). 
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50  second  group  velocities 


20  second  group  velocities 


Figure  7.  Map  of  group  velocities  at  50  seconds  (top)  and  20  seconds  (bottom). 


13 


Figure  8  shows  phase  velocity  as  a  function  of  frequency  for  all  models,  and  Figure  9  shows 
a  map  of  phase  velocity  at  50  and  20  seconds. 
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Figure  8  Average  phase  velocity  plotted  vs.  frequency  and  ±  one  standard  deviation.  Top 
figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and  oceanic  models 
(blue). 
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50  second  phase  velocities 


20  second  phase  velocities 


Figure  9.  Map  of  phase  velocities  at  50  seconds  (top)  and  20  seconds  (bottom). 
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Figure  10  shows  the  attenuation  coefficient,  gamma,  as  a  function  of  frequency  for  all  models, 
and  Figure  1 1  shows  a  map  of  gamma  at  50  and  20  seconds.  The  attenuation  coefficient  is 
derived  from  generic  models,  not  through  inversion  of  attenuation  data,  so  should  be  regarded  as 
approximate. 
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Figure  10  Average  attenuation  coefficient  plotted  vs.  frequency  and  ±  one  standard  deviation. 
Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and  oceanic 
models  (blue). 
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50  second  gamma 


20  second  gamma 


Figure  11.  Map  of  gamma  at  50  seconds  (top)  and  20  seconds  (bottom). 
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Figure  12  shows  the  source  region  excitation  function,  SI,  as  a  function  of  frequency  for  all 
models,  and  Figure  13  shows  a  map  of  SI  at  50  and  20  seconds. 
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Figure  12  Average  source  amplitude,  SI,  plotted  vs.  frequency  and  ±  one  standard  deviation. 
Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and  oceanic 
models  (blue). 


18 


50  second  si  source  amplitude 


20  second  si  source  amplitude 


Figure  13.  Map  of  S 1  at  50  seconds  (top)  and  20  seconds  (bottom). 


19 


Figure  14  shows  the  receiver  region  amplitude  function,  S2,  as  a  function  of  frequency  for  all 
models,  and  Figure  1 5  shows  a  map  of  S2  at  50  and  20  seconds. 
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Figure  14  Average  source  amplitude,  S2,  plotted  vs.  frequency  and  ±  one  standard  deviation. 
Top  figure  is  for  all  earth  models.  Bottom  figure  shows  continental  models  (red)  and  oceanic 
models  (blue). 
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50  second  s2  receiver  amplitude 


Figure  15.  Map  of  S2  at  50  seconds  (top)  and  20  seconds  (bottom). 
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4  OPTIMIZED  SURFACE  WAVE  AMPLITUDE  MEASUREMENTS 


Surface  wave  measurements  traditionally  are  accomplished  by  measuring  a  time  domain 
amplitude  at  a  period  near  20  seconds  and  then  calculating  a  surface  wave  magnitude  Ms.  This 
procedure  is  problematic  at  regional  distances  because  the  surface  wave  is  not  well  dispersed  and 
a  distinct  20-second  arrival  may  not  be  present.  It  is  possible  to  measure  time  domain  amplitudes 
at  higher  frequencies  with  corrections  (e  g.  Marshall  and  Basham,  1972),  however  measurements 
may  be  inaccurate  due  to  differences  in  dispersion  caused  by  differences  in  earth  structure. 
Stevens  and  McLaughlin  (2001)  suggested  as  an  alternative  replacing  time  domain 
measurements  with  a  path  corrected  spectral  magnitude. 

4.1  Path  Corrected  Spectral  Magnitudes 

The  path  corrected  spectral  magnitude,  logM0,  is  calculated  by  dividing  the  observed  surface 
wave  spectrum  by  the  Green’s  function  for  an  explosion  of  unit  moment  and  taking  the  logarithm 
of  this  ratio,  averaged  over  any  desired  frequency  band.  The  path  corrected  spectral  magnitude  is 
defined  as  the  logarithm  of: 


M0  = 


U(a>,r,0)/ 


St (a,K )S2(co) exp [- v  (co )r] A 


Jats\n(r  I  at) 


(4.1) 


where  U  is  the  observed  surface  wave  spectrum,  and  as  described  earlier  depends  on  the 
source  region  elastic  structure  and  the  explosion  source  depth,  S2  depends  on  the  receiver  region 

elastic  structure,  yp  is  the  attenuation  coefficient  that  depends  on  the  attenuation  integrated  over 
the  path  between  the  source  and  receiver,  r  is  the  source  to  receiver  distance  and  ae  the  radius  of 
the  earth.  All  of  the  functions  in  equation  4. 1  are  easily  derived  from  plane-layered  earth  models, 
and  allow  the  measurement  to  be  regionalized  to  account  for  differences  in  earth  structure  at  the 
source  and  receiver,  and  due  to  attenuation  along  the  path. 

The  advantages  of  using  logM0  instead  of  the  traditional  surface  wave  magnitude  A/s  are  that 
logATo  is  insensitive  to  dispersion,  independent  of  distance,  works  well  at  regional  distances,  and 
can  be  regionalized.  Regionalized  path  corrected  spectral  magnitudes  incorporate  geographic 
variations  in  source  excitation  and  attenuation.  Furthermore,  as  discussed  below,  it  can  in 
principle  be  measured  over  different  frequency  bands  to  optimize  the  signal-to-noise  ratio.  Mt 
and  logMo  share  some  limitations:  spectra  from  earthquakes  vary  due  to  source  mechanism  and 
depth,  and  errors  can  occur  if  the  measurement  is  made  in  a  spectral  dip  or  at  high  frequencies 
for  a  deep  event  (Figure  16).  Azimuthal  variations  in  amplitude  caused  by  focal  mechanism  also 
affect  the  amplitudes  of  both  logA70  and  Ms.  logMy  can  also  be  corrected  for  structural 
heterogeneity  using  the  amplitude  corrections  described  earlier. 
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Figure  16.  Path  corrected  spectra  for  an  explosion  and  for  earthquakes  calculated  for  several 
depths.  The  path  corrected  explosion  spectrum  is  flat  over  the  entire  frequency  band  (for  perfect 
data  and  path  correction),  while  the  path  corrected  earthquake  spectrum  is  flattened,  but  has 
some  variation  due  to  source  mechanism  and  source  depth. 


The  test  cases  discussed  by  Stevens  and  McLaughlin  (2001)  used  a  frequency  band  of  0.02-0.05 
Hz  (50-20  second  period)  to  estimate  the  spectral  magnitudes.  They  estimated  that  on  average, 
the  time  domain  and  spectral  magnitudes  are  related  as  logMo=M,+ 1 1 .75.  Most  of  the  waveforms 
in  that  work  were  recorded  at  distances  exceeding  8°  Due  to  the  relatively  flat  path  corrected 
spectra  over  the  0.02-0.05  Hz  band  for  most  data,  this  choice  worked  quite  well.  The  authors 
noted,  however,  that  higher  frequencies  might  be  required  for  shorter  paths.  An  important 
observation  was  that  the  logM0  residuals  are  independent  of  distance,  despite  the  simple  O 
models  used  in  the  earth  structures  (Figure  17). 
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Figure  17.  Path  corrected  spectral  magnitude  (Log  Mo)  residuals  plotted  vs.  distance  (from 
Stevens  and  McLaughlin,  2001).  Log  M0  is  nearly  independent  of  distance. 
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In  the  following  test  case,  we  look  at  the  utility  of  higher  frequencies  for  estimating  spectral 
magnitudes  of  smaller  events,  recorded  at  smaller  distances.  The  purpose  is  to  optimize  the 
spectral  magnitude  estimates,  to  test  their  distance  and  frequency  independence,  and  to  identify 
any  measurement  problems  or  pitfalls.  For  large  amplitude  signals  we  can  expect  the  lower 
frequencies  to  be  better  in  general,  particularly  at  larger  distances  due  to  greater  attenuation  at 
higher  frequencies.  Our  hypothesis  at  the  initiation  of  the  study  was  that  using  higher  frequencies 
for  measuring  spectral  magnitudes  at  shorter  distances  would  optimize  signal  to  noise  ratio  and 
therefore  be  better  for  measuring  surface  waves  at  regional  distances,  however  as  discussed 
below  this  is  only  true  to  a  limited  extent  and  may  come  at  the  cost  of  degraded  discrimination. 

4.2  Surface  Wave  Measurements  Using  Events  On  and  Near  the  Lop  Nor  Test  Site 

To  optimize  the  measurement  procedures  and  examine  the  performance  of  the  path  corrected 
spectral  magnitude  at  regional  distances,  we  used  584  spectra  from  76  earthquakes  and  1 1 
explosions  in  the  Lop  Nor  area  (Figure  18).  Approximately  1 1%  of  the  spectra  used  for  the 
logMo  estimates  originate  from  records  at  source-station  distances  of  5°  or  less,  and  another  1 1% 
at  distances  of  30°  or  greater.  The  bulk  of  the  data  comes  from  intermediate  distances. 


50‘ 
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Figure  18.  Maps  showing  the  Lop  Nor  area  (rectangle),  stations  (triangles),  and  earthquake 
(circles)  and  explosion  (crosses)  epicenters. 

Figure  19  shows  examples  of  explosion  and  earthquake  path  corrected  spectra  from  Lop  Nor  at 
various  distances.  The  explosion  spectra  are  flattened  as  expected,  while  there  is  more  variation 
of  the  earthquake  spectra  with  frequency.  This  happens  because  1 )  the  spectra  are  corrected  by 
an  explosion  Green’s  function  that  flattens  earthquake  spectra  imperfectly,  and  2)  the  earthquake 
spectra  have  frequency  variations  caused  by  source  mechanism  and  depth.  Most  of  the  spectra 
shown  in  Figure  19  are  contaminated  by  noise  at  the  lowest  frequencies. 
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Figure  19.  Examples  of  path  corrected  spectra  used  in  this  work:  (a)  Lop  Nor  explosions  recorded 
at  distances  of  2°,  7°  and  67°  (left);  (b)  Lop  Nor  earthquakes  recorded  at  distances  of  0.4°,  22°  and 
65°  (right).  See  examples  of  station  logM0  estimates  in  Table  1 .  S/N  is  good  at  all  but  the  highest 
frequencies.  Low  frequency  noise  dominates  over  the  surface  wave  signal  below  about  0.02  Hz. 


We  calculated  individual  spectral  magnitudes  (i.e.,  several  station  magnitudes  per  event)  over  all 
possible  frequency  bands  between  0.02  Hz  (50  s)  and  0. 15  Hz  (~7s),  with  bandwidths  of  0.03 
Hz,  0.04  Hz,  etc.,  up  to  0. 1 3  Hz  for  the  0.02-0. 1 5  Hz  band.  This  procedure  provided  1 53  bands  to 
examine  from  each  spectrum.  Four  different  methods  were  used  as  follows  to  identify  the  most 
robust  spectral  magnitude  estimate 

1 .  Calculating  a  “simple”  mean  of  the  logarithms  of  all  path  corrected  amplitude 
measurements  made  in  a  given  frequency  band.  This  is  comparable  to  Stevens  and 
McLaughlin’s  (2001)  estimates  in  the  0.02-0.05  Hz  band. 

2.  Iteratively  calculating  a  “robust”  mean,  by  rejecting  outliers  outside  two  standard 
deviations  from  the  mean  calculated  at  each  step.  The  procedure  ends  when  all 
measurements  remain  within  two  standard  deviations  or  when  half  of  the  amplitude 
measurements  in  a  frequency  band  are  rejected,  whichever  occurs  first.  Thus  the  spectral 
magnitude  estimates  are  much  less  affected  adversely  by  the  tendency  of  some  spectra  to 
sharply  vary  in  amplitude  over  some  frequencies,  with  most  outliers  marking  anomalously 
low  amplitudes  (see  Figure  19  above).  Figure  20  compares  the  individual  (station)  log/V/o 
estimates  from  (1)  and  (2).  Standard  deviations  from  the  robust-mean  method  are 
predictably  lower  than  those  in  the  simple-mean  method,  as  the  insets  in  Figure  20  show. 

3.  Calculating  logMo  at  the  center  frequency  of  a  least-squares  straight  line  fit  to  the 
spectrum  over  a  given  frequency  band. 

4.  Same  as  (3),  but  the  straight  line  is  “robust”,  minimizing  the  absolute  deviations  of  the 
logarithms  of  the  amplitudes  from  the  line. 
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Figure  20.  Comparison  of  station  spectral 
magnitudes  calculated  with  two  different 
methods.  Insets  show  standard  deviations  as 
indicated.  See  text  for  details. 


Figure  21.  Histograms  of  standard  deviations  of 
the  mean  station  logMi  estimated  in  three 
frequency  bands  (shown  on  right),  for  small  and 
larger  distances  (shown  on  top). 


The  above  estimates  were  compared  in  order  to  select  the  most  suitable  frequency  bands, 
possibly  varying  with  distance.  Ideally,  the  corrected  spectra  would  be  flat  over  an  extended 
band  of  frequencies.  Flatness  is  particularly  expected  for  explosions,  as  supported  (within  limits) 
by  the  explosion  examples  in  Figure  19  above.  The  magnitude  spectra  estimated  over  any 
reasonable  band  would  be  then  consistent.  In  reality,  completely  flat  spectra  over  extended 
frequency  bands  are  rare,  so  we  need  to  choose  bands  small  enough  not  to  include  too  many 
variable  features  of  the  spectra,  yet  large  enough  not  to  reflect  only  local,  possibly  spurious, 
characteristics. 

In  view  of  the  above,  the  two  main  desirable  properties  of  a  spectrum  over  a  given  frequency 
band  are  small  standard  deviations  and  flatness.  For  this  reason,  in  our  search  for  optimum 
frequency  bands  we  used  two  criteria.  First,  small  standard  deviations  from  (1)  above  represent 
one  measure  of  the  suitability  of  a  frequency  band.  Figure  21  indicates  that  for  small  distances, 
the  0.08-0. 1 1  Hz  frequency  band  may  be  preferable  (the  largest  number  of  small  standard 
deviations)  to  either  0.02-0.05  Hz,  or  0. 12-0. 1 5  Hz.  Larger  distances  do  not  present  a  clear 
picture,  but  it  is  still  evident  that  relatively  more  small  standard  deviations  are  found  in  the  0.02- 
0.05  Hz  frequency  band,  compared  with  the  higher  frequencies.  We  note  that  at  this  stage  we  do 
not  use  the  standard  deviations  from  (2),  since  they  are  designed  to  greatly  diminish  the  presence 
of  outliers  and  are  thus  not  representative  enough  of  the  quality  of  the  estimates  in  the  different 
frequency  bands.  However,  once  a  suitable  frequency  band  is  chosen,  the  robust  mean  is  the 
most  reliable  estimate  of  logMo.  Spectral  flatness  as  measured  with  the  slopes  of  the  “robust” 
lines  in  (4)  above  provides  a  second  measure  of  the  quality  of  frequency  band;  the  smaller  the 
slope,  the  flatter  the  spectra.  Table  1  shows  examples  of  selected  estimates,  over  one  specific 
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frequency  band  out  of  1 53  (0.05-0. 1  Hz),  for  the  explosion  and  earthquake  spectra  shown  in 
Figure  19.  Smaller  slopes  (flatter  spectra)  are  evident  for  explosions  compared  with  earthquakes. 
On  the  other  hand,  increasing  absolute  values  of  slopes  and  standard  deviations  are  seen  for 
earthquakes  with  increasing  distance.  This  is  to  be  expected,  as  the  relatively  high-frequency 
band  in  the  example  is  less  suitable  as  distance  increases. 


Table  1.  Station  LogMO  estimates  in  0.05-0. 10  Hz  from  the  spectra  in  Figure  19. 


Event  Type 

ID.  station 

Distance 
,  degrees 

Station  logM{) 
(simple) 

Station  logMo 
(robust) 

Station 

slope/logM0 

mh 

Explosions 

21450528. W 
MQ 

2.2 

14.31+0.10 

14.33+0.08 

+1.02/14.31 

4.5 

21450535. M 
AK 

7.1 

15.60+0.15 

15.67+0.05 

+0.41/15.64 

5.8 

21450534.ES 

DC 

66.8 

14.95+0.22 

15.02+0.09 

+1.62/14.97 

5.4 

Earthquakes 

21456615. W 
MQ 

0.4 

13.92+0.14 

14.01+0.05 

+3.64/13.93 

3.2 

21456712. AR 
U 

22.2 

14.44+0.25 

14.44+0.24 

-8.70/14.48 

3.8 

21457058.IL 

AR 

65.3 

15.55+0.27 

15.45+0.10 

-11.18/15.59 

4.3 

14  15  16 

LogMO  (0.02-0.05  Hz) 


Figure  22.  Comparison  of  station  spectral  magnitudes  in  different  frequency  bands:  (a)  adjacent 
bands,  all  distances;  (b)  overlapping  bands,  all  distances;  (c)  overlapping  bands,  distances  <5°. 


Next,  we  examined  the  consistency  of  spectral  magnitude  estimates  in  different  frequency  bands. 
Figure  22  shows  examples  of  such  estimates  in  several  frequency  bands  (marked  along  the  plot 
axes).  These  results  indicate  that  although  measurements  are  generally  consistent  when  measured 
in  different  frequency  bands,  some  individual  measurements  do  change  significantly.  Also,  there 
is  a  tendency  for  measurements  to  be  smaller  at  higher  frequencies  (points  lie  slightly  to  the  right 
of  the  lines  in  Figure  22).  These  results  indicate  that  spectral  magnitudes  can  be  measured  in 
different  frequency  bands,  but  with  some  caution  and  attention  to  spectral  shape  variations. 
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We  examined  which  frequency  bands  performed  best  for  discrimination  between  small 
earthquakes  and  explosions.  That  is,  we  evaluated  whether  any  set  of  variable  frequencies 
discriminate  better  than  a  single  frequency  band  applied  at  all  distances.  Figure  23  (left)  shows 
log/Wofl/b  plots  using  a  set  of  variable  frequencies  (0.02-0.05  Hz  for  distances  exceeding  25°, 
0.06-0.09  Hz  for  10°  to  25°,  and  0.08-0. 1 1  Hz  for  distances  below  10°).  In  a  large  statistical 
analysis  of  measurements  made  in  different  frequency  bands  at  different  distances,  these 
frequency  bands  were  found  to  be  the  optimum  for  measurement,  based  on  best  signal  to  noise 
and  flat  spectral  shape.  Figure  23  (right)  shows  logMiWb  plots  using  the  0.03-0.07  Hz  frequency 
band  for  all  distances. 

The  plot  on  the  left,  where  higher  frequencies  are  used  at  small  distances  (and  hence  for  the 
smallest  earthquakes)  appears  to  discriminate  less  well  for  small  events  than  when  a  constant 
0.03-0.07  Hz  frequency  band  is  used  for  the  magnitude  measurements.  We  examined  the 
logMo:wb  ratio  for  a  number  of  frequency  bands  and  established  that  the  0.03-0.07  Hz  interval 
performs  best  for  discriminating  between  earthquakes  and  explosions  for  the  Lop  Nor  data  set. 
This  effect  might  be  expected  from  the  path  corrected  spectra  shown  in  Figure  16,  since  the 
earthquake  spectra  exhibit  a  slow  decline  in  amplitude  with  frequency.  This  experiment 
demonstrates  a  potentially  important  tradeoff  -  for  small  events  higher  frequencies  may  provide 
the  best  signal  to  noise  performance,  but  have  degraded  discrimination  capability. 


Figure  23.  Log /V/0:/wb  plots  showing  event  spectral  magnitudes  for  earthquakes  (O)  and 
explosions  (X)  in  Lop  Nor.  Station  spectral  magnitudes  were  calculated  using  frequencies 
increasing  with  decreasing  distance  (left)  and  the  0.03-0.07  Hz  frequency  band  for  all  distances 
(right).  Bold  lines  mark  the  empirical  discrimination  relationship  of  Stevens  and  McLaughlin 
(2001).  For  small  events  higher  frequencies  may  provide  the  best  signal  to  noise  performance, 
but  at  a  cost  of  degraded  discrimination  capability. 


4.3  Path  Corrected  Noise  Estimates  at  the  Lop  Nor  Test  Site 

While  a  path  corrected  spectral  magnitude  can  be  measured  over  any  frequency  band,  it  is 
subject  to  the  following  constraints: 
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1 .  Earthquake  spectra  decrease  at  high  frequencies,  depending  on  depth  (see  Figure  16),  so 
the  high  end  of  the  frequency  band  should  be  low  enough  that  discrimination  is  not 
degraded  by  this  effect. 

2.  At  high  frequencies,  attenuation  is  greater  and  the  dispersion  more  variable,  so  the  path 
correction  is  likely  to  be  less  well  determined  and  the  signal  may  be  smaller  than  at  lower 
frequencies,  particularly  at  longer  distances. 

3.  Noise  increases  at  low  frequency,  so  the  low  end  of  the  frequency  band  should  be  high 
enough  to  be  above  the  noise  level. 

In  order  to  quantify  these  effects,  we  created  some  “path  corrected  noise  spectra”  using  path 
corrections  appropriate  for  the  Lop  Nor  test  site.  These  are  simply  noise  spectra  measured  at  the 
station  that  have  been  divided  by  an  explosion  Green’s  function  in  the  same  manner  as  would  be 
done  for  a  signal.  Since  the  signal  spectra  are  approximately  flat  over  most  of  the  frequency 
band,  the  path  corrected  noise  spectra  are  a  measure  of  the  minimum  path  corrected  signal  that 
could  be  measured  at  each  station. 


Path  Com**  Mo.v*  Spadra  at  WMO 


Path  CorracMd  Nom*  Spaeha  at  HtA 


Figure  24.  Average  “path  corrected  noise  measurement”  and  ±one  standard  deviation  curves  for 
13  time  segments  at  WMQ  (left)  and  for  54  time  segments  at  H1A  (right).  The  path  corrected 
noise  measurement  gives  the  minimum  path  corrected  signal  level  that  could  be  observed  at  each 
frequency.  Recall  that  the  predicted  path  corrected  explosion  spectrum  is  flat  over  this  frequency 
band,  and  the  predicted  path  corrected  earthquake  spectra  decrease  with  increasing  frequency 
(see  Figure  16). 

Figure  24  shows  the  path  corrected  noise  measurements  for  two  stations:  WMQ,  located  an 
average  of  2  degrees  from  each  event,  and  HIA,  located  an  average  of  27  degrees  from  each 
event.  At  HIA,  the  noise  spectrum  has  minima  at  0.04  and  0. 1  Hz,  either  side  of  the  0.06  Hz 
primary  microseism  peak.  Noise  amplitudes  increase  substantially  below  about  0.04  Hz.  Noise 
levels  at  WMQ  are  similar,  but  are  flatter  from  0.04  to  0.2  Hz  with  a  less  prominent  microseism 
peak  and  lower  (path  corrected)  noise  at  the  highest  frequencies  because  of  the  close  distance.  As 
Figure  24  shows,  there  is  also  considerable  variability  in  the  noise  level,  so  that  although  the 
average  noise  level  is  about  13.5  at  WMQ  and  14.5  at  HIA  (corresponding  approximately  to  M, 
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1.7  and  2.7,  respectively),  the  minimum  spectral  magnitude  that  could  be  measured  may  be 
substantially  lower  or  higher  depending  on  the  noise  level  at  the  time  the  measurement  is  made. 

These  results  suggest  that  in  general  it  is  best  to  measure  surface  wave  spectra  at  frequencies 
above  0.03  Hz  and  below  0. 1  Hz.  Although  noise  levels  remain  fairly  low  in  the  0. 1-0.2  Hz 
frequency  band  at  the  closest  stations,  factors  #1  and  #2  mentioned  above  make  measurement  in 
this  band  risky.  Based  on  a  number  of  empirical  tests,  we  recommend  measuring  surface  waves 
at  frequencies  between  0.03  and  0.08  Hz  for  this  region.  Further  research  is  required  to  determine 
the  optimum  frequency  band  in  other  regions  where  there  are  structures  that  cause  large 
amplitude  changes,  such  as  along  paths  that  cross  oceanic  regions. 

4.4  Comparison  of  Spectral  and  Time  Domain  Magnitudes 

Russell  (2004)  proposed  a  new  type  of  surface  wave  magnitude  M,(b)  which  differs  from  a 
traditional  20  second  magnitude  in  that  it  uses  a  Butterworth  filter  to  measure  a  time  domain 
amplitude  in  a  narrow  band  around  any  desired  frequency,  and  then  applies  a  correction  for  the 
source  function  similar  to  the  explosion  source  function  used  in  the  path  corrected  spectral 
magnitude  described  above.  The  main  purpose  of  Ms(b)  is  to  allow  surface  waves  to  be  measured 
at  regional  distances  at  higher  frequencies.  The  magnitude  is  defined  by: 


Ms(b)  =  log  {Ab )  +  ~  log  (sin  A)  +  0.003 1 


20^ 

(  20^ 

A -0.66  log  — 

T  j 

\T  ) 

-log(/c)-0.43 


(4.2) 


where  Ab  is  the  filtered  amplitude,  T  is  the  measured  period,  and  fc  is  the  Butterworth  filter 
width.  It  is  instructive  to  compare  the  terms  in  the  Russell  magnitude  with  the  Rezapour  and 
Pearce  (1998)  Ms  and  the  path  corrected  spectral  magnitude  logAT0  described  above.  This  is 
shown  in  Table  2. 


Table  2.  Comparison  of  time  domain  and  spectral  magnitude  measurement  and  correction  terms. 
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Note  that  each  magnitude  makes  a  slightly  different  set  of  corrections.  logM0  corrects  for  both 
source  and  receiver  structure  based  on  earth  models  at  those  locations  and  using  an  explosion 
Green’s  function  at  the  source.  Similarly,  Ms(b>  applies  a  source  correction  based  on  typical 
explosion  source  excitation.  The  Airy  phase  dispersion  correction  accounts  for  superposition  of 
waves  with  similar  group  velocities,  and  is  needed  only  in  the  time  domain.  Ms(b)  uses  a 
Butterworth  filter  that  is  sufficiently  narrow  to  avoid  this  problem.  The  filter  correction  corrects 


30 


for  the  width  of  the  Butterworth  filter.  The  normalization  for  the  two  Ms  measurements  is  chosen 
to  make  it  consistent  with  historical  Ms  magnitudes  at  a  chosen  distance  range.  logMo  has  natural 
units  of  log  moment  and  is  not  otherwise  normalized,  however  Stevens  and  McLaughlin  (2001 ) 
showed  that  logM0  -1 1 .75  is  consistent  with  the  Rezapour  and  Pearce  Ms.  Attenuation  for  Ms  is 
an  empirical  correction  based  on  a  very  large  number  of  20  second  measurements.  Ms(b)  similarly 
uses  an  empirical  attenuation  correction,  but  also  includes  an  empirical  correction  for  the  change 
in  attenuation  with  frequency.  logM0  uses  attenuation  calculated  from  earth  (velocity,  density 
and  Q)  models  along  a  source  to  receiver  path. 

A  path  corrected  time  domain  magnitude  can  be  derived  by  combining  the  path  corrected 
spectral  magnitude  with  M,(b),  using  the  source  and  path  corrections  from  earth  models  to  replace 
the  empirical  average  corrections.  We  define  the  path  corrected  time  domain  magnitude  Ms{bP)  as: 

Ms(bP)  =  *o§  M  +  \  log  (sin  A)  +  M  "  ,08  )-  log(52 )  -  log  (fc )  +  Cbp  (4.3) 

where  Cbp  is  a  constant  chosen  to  make  Ms(bP)  consistent  with  historical  magnitudes.  Although 
equation  4.3  may  appear  more  complicated  than  equation  4.2,  the  functions  Si,  S2,  and  yp  are 
easily  tabulated  and  stored  in  files,  and  a  computer  can  quickly  calculate  them  for  any  path  based 
on  a  simple  lookup  table.  There  is  substantial  regional  variation  in  these  quantities  that  should  be 
removed  to  ensure  consistent  measurements  (examples  of  Si,  S2,  and  y  for  continental  and 
oceanic  structures  are  shown  in  Stevens  and  McLaughlin,  1996).  Another  advantage  of  this 
approach  is  that  it  can  ensure  that  fc,  which  must  be  less  than  a  minimum  value  calculated  from 
the  group  velocity,  is  always  set  appropriately. 

Bonner  et  al.  (2004)  showed  that  M,(b>  gave  consistent  results  as  a  function  of  distance,  and  also 
gave  good  discrimination  results  when  the  maximum  value  of  M«(b)  over  the  entire  frequency 
band  was  used.  This  result  could  be  expected  from  Figure  16,  showing  the  frequency  variation 
of  earthquake  spectra  as  a  function  of  depth,  and  from  Figure  19,  showing  the  variability  of 
measured  path  corrected  earthquake  spectra.  By  taking  the  maximum  value  of  Ms(b)  instead  of 
measuring  all  surface  waves  at  the  same  frequency,  Ms(b)  can  move  to  lower  frequencies  for 
earthquakes,  avoiding  the  spectral  reduction  at  higher  frequencies,  and  can  avoid  spectral  dips 
that  commonly  occur  in  real  data.  This  may,  therefore,  be  a  better  approach  to  measuring  surface 
waves  for  discrimination  than  the  distance  dependent  or  fixed  frequency  band  spectral  methods 
described  above.  However,  some  caution  is  necessary  since  as  shown  in  Figure  24,  noise 
increases  rapidly  at  lower  frequencies,  so  the  maximum  value  could  easily  correspond  to  a  noise 
measurement  even  when  a  surface  wave  is  clearly  visible  at  higher  frequencies. 
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5 


BACKAZIMUTH  ESTIMATION  RELIABILITY  USING  SURFACE 
WAVE  POLARIZATION 


As  part  of  its  responsibilities  under  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT),  the 
International  Data  Center  (1DC)  automatically  processes  seismic  data  recorded  at  International 
Monitoring  System  (IMS)  stations,  identifying  surface  waves  and  measuring  surface  wave 
amplitude  and  Ms  (Stevens  and  McLaughlin,  2001).  The  Ms  measurement  is  important  because 
of  the  earthquake/explosion  discrimination  capability  of  the  M9:mb  discriminant  (e  g.  Marshall 
and  Basham,  1972;  Stevens  and  Day,  1985).  Three-component  backazimuth  estimates  derived 
using  the  current,  spectral  method  are  recorded  but  are  too  inaccurate  to  use  for  association.  That 
method  is  based  on  an  algorithm  originally  proposed  (Smart,  1978)  as  part  of  a  surface  wave 
detector,  rather  than  as  a  means  of  estimating  the  backazimuth.  Selby  (2001 )  suggested  that 
another  detection  technique  (Chael,  1997),  based  on  cross-correlation  of  the  vertical  and  Hilbert 
transformed  radial  component  Rayleigh  wave  records,  could  be  used  for  backazimuth  estimation 
and  would  improve  IDC  processing. 

To  compare  the  performance  of  the  algorithms,  we  incorporated  the  Chael/Selby  (CS)  algorithm 
into  the  automatic  surface  wave  processing  software  that  utilizes  the  current  IDC  method  (CM). 
Under  current  IDC  procedures,  for  an  event  location  and  time  based  on  body  wave  arrivals,  a 
surface  wave  is  considered  detected  if  it  passes  a  dispersion  test  in  the  group  velocity  window 
predicted  for  the  event.  We  test  both  algorithms  on  2,599  records  that  passed  the  dispersion  test 
and  2,363  generally  poorer  S/N  records  that  failed,  all  from  events  that  had  at  least  four  surface 
wave  detections.  We  assess  the  performance  of  each  algorithm  for  different  passbands  and  group 
velocity  windows. 

5.1  Differences  Between  the  Azimuth  Estimation  Algorithms 

Both  the  CM  and  CS  algorithms  find  the  backazimuth  that  best  matches  the  expectation  that  the 
horizontal  and  vertical  components  of  the  Rayleigh  wave  records  are  similar,  but  90°  out  of 
phase.  The  biggest  difference  between  algorithms  is  that  CM  uses  Love  as  well  as  Rayleigh 
waves  to  estimate  the  backazimuth.  Consequently,  the  CS  algorithm  should  perform  optimally 
when  the  time  window  encompasses  just  the  Rayleigh  wave,  as  a  longer  window  will  only  add 
noise,  while  the  CM  should  perform  optimally  with  a  longer  time  window  that  encompasses  both 
the  Rayleigh  wave  and  the  higher  group  velocity  Love  wave. 

5.1.1  The  C  urrent  Method 

The  CM  finds  values  for  rn  and  In,  the  complex  Fourier  coefficients  of  the  Rayleigh  and  Love 
wave  vertical  component  displacements,  <|>,  the  backazimuth,  and  e„,  the  Rayleigh  wave 
ellipticity,  that  minimize  the  squared  distance  between  the  data  and  model  (Smart,  1978;  Smart 
and  Sproules,  1981).  The  ellipticity  may  be  assumed  known,  which  reduces  by  one  the  degrees 
of  freedom.  Subscripts  refer  to  frequency.  The  function  minimized  is 

Z  \zn~rn\2 +\yn-{isnrna-lnfif +\xn-(ienrn/3  +  lna)[  (5.1) 
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where  xn,yn,  and  z„  are  the  complex  Fourier  coefficients  of  the  east,  north,  and  vertical 
components  of  the  seismic  record  at  the  nth  frequency,  and  a  and  p  are  cos(4>)  and  sin(<|)) 
respectively.  Two  important  elements  of  this  function  are  that  the  radial  component  of  the 
Rayleigh  wave  equals  ie„rn,  that  is,  the  vertical  component  advanced  by  90°  and  scaled  by  the 
ellipticity,  and  that  the  Love  wave  is  independent  of  the  Rayleigh  wave. 

5. 1 .2  The  Chael/Selby  Algorithm  and  its  Implementation 

The  CS  algorithm  (Chael,  1997;  Selby,  2001)  finds  the  backazimuth  for  which  the  correlation 
coefficient  of  the  vertical  and  Hilbert  transformed  radial  is  a  maximum  (Equation  5.2). 


(5.2) 


where. 


S;*=X*,(rK(r)  (53) 

r=l 

The  implementation  uses  the  cross-correlation  of  the  Hilbert  transformed  vertical  with  the  radial, 
for  one  degree  increments  of  backazimuth.  While  Selby  uses  the  maximum  positive  cross¬ 
correlation  as  the  backazimuth,  Chael  (1997)  uses  the  central  azimuth,  determined  by  the  circular 
mean  (e  g.  Fisher,  1993).  We  test  the  algorithm  using  both  measures. 

For  synthetic  Rayleigh  waves  with  azimuthally  evenly  distributed  random  noise,  the  circular 
mean  provides  more  accurate  estimates  than  the  maximum  cross-correlation.  The  lower  the  S/N, 
the  greater  the  advantage  of  the  circular  mean. 

For  noise  free  data  (i.e.  synthetics),  equation  5.2  will  return  a  negative  or  positive  constant,  as  the 
numerator  changes  in  synch  with  the  second  term  of  the  denominator  (the  first  denominator  term 
is  constant)  as  the  algorithm  steps  through  backazimuth.  Our  implementation  of  the  algorithm 
therefore  provides  two  methods  of  estimating  the  similarity  of  the  components.  One  uses  the 
correlation  coefficient  (equation  5.2).  The  second  method  is  intended  to  avoid  the  problem 
described  above  when  data  are  noise  free  by  normalizing  by  S::  alone  (equation  5.4). 


5.2  Performance  of  the  Algorithms  in  4  Frequency  Bands 

The  first  set  of  tests  used  a  2.5  to  5.0  km/sec  group  velocity  window,  which  should  encompass 
both  Love  and  Rayleigh  waves,  and  compared  backazimuth  estimates  in  3  relatively  narrow, 
overlapping  frequency  bands,  0.02-0.04  Hz,  0.03-0.05,  0.04-0.06  Hz,  plus  one  frequency  band 
covering  the  entire  spectrum  from  0.02  to  0.06  Hz. 

Figure  25  shows  histograms  of  backazimuth  residuals  for  the  CM  and  one  implementation  of  the 
CS  algorithm  for  the  middle  frequency  band.  Table  3  presents  the  errors  associated  with  each 
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method  The  CS  algorithm  performs  better  in  two  ways.  First,  the  CM  applied  as  it  is  currently 
used  in  the  automatic  processing,  has  errors  of  180°  for  a  significant  number  of  events.  Second, 
the  histograms  reveal  many  more  outliers  for  the  CM  than  for  the  CS  algorithm. 
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Figure  25.  Backazimuth  residuals  for  0.03-0.05  Hz  using  the  current  method  (left)  and  the 
Chael/  Selby  algorithm  using  the  normalization  of  equation  5.2  and  peak  cross-correlation 
(right). 


Surprisingly,  the  best  implementation  of  the  CS  algorithm  is  the  combination  of  the  peak 
position  of  the  cross-correlation  and  normalization  of  the  correlation  by  yfsls^  Why  the 

maximum  cross-correlation  provides  a  more  accurate  estimate  of  backazimuth  than  the  circular 
mean  bears  further  investigation.  One  possibility  is  that  the  real  noise  is  azimuthally  unevenly 
distributed,  consistent  with  observations  of  Selby  (2001 ),  and  so  biases  the  circular  mean 
measurement  toward  the  direction  of  the  predominant  noise  source.  Such  skewed  noise  would 
not  so  strongly  affect  the  maximum  cross-correlation  position. 

We  use  the  scaled  median  absolute  deviation  (SMAD),  the  one-norm  measure  of  the  central 
tendency,  to  provide  a  measure  of  the  spread  of  data  about  the  central  value  that  is  less  biased  by 
outliers  than  the  STD  in  heavy-tailed  data  such  as  these. 

The  CM  performs  best  at  higher  frequency,  while  the  CS  algorithm  performs  best  at  lower 
frequency.  The  CM  performs  almost  as  well  as  the  CS  algorithm  at  the  highest  frequencies 
tested.  This  is  likely  because  the  Love  waves  are  relatively  larger  at  longer  periods  than  the 
Rayleigh  waves,  and  contrary  to  expectations,  the  CM  performs  poorly  when  the  Love  waves  are 
very  large.  This  is  due  to  the  180  degree  ambiguity  in  Love  wave  polarization,  which  causes 
large  Love  wave  amplitudes  to  increase  the  likelihood  of  a  180  degree  error  in  backazimuth. 
Larger  backazimuth  residuals  at  higher  frequency  for  the  CS  method  are  due  to  lower  S/N. 

For  the  CS  algorithm,  the  backazimuth  residuals  for  the  broader  frequency  band  are  comparable 
to  those  in  the  0.03-0.05  Hz  passband,  larger  than  those  of  the  lowest  frequency  passband,  and 
smaller  than  those  of  the  highest  frequency  passband.  The  parts  of  the  frequency  band  where  the 
backazimuth  either  varies  or  is  less  well  resolved  diminish  resolution  of  the  broadband  estimate 


34 


compared  to  narrower  passband  estimates,  so  the  broader  frequency  band  does  not  provide  an 
advantage  For  the  CM,  the  performance  is  poorer  for  the  broadband  data  than  for  each  of  the 
narrower  passbands  used.  The  rest  of  the  analyses  are  performed  for  just  the  three  narrow 
frequency  bands. 


Table  3.  Scaled  median  absolute  deviations  and  standard  deviations  (in  parentheses)  of  the 
backazimuth  residuals  in  4  frequency  bands  from  the  algorithms  tested.  Four  variations  of  the 
implementation  of  the  CS  algorithm  are  tested.  CSi  and  CS2  use  the  peak  amplitude  of  the  cross¬ 
correlation  while  CS3  and  CS4  use  the  circular  mean.  CSi  and  CS3  are  normalized  by  \jS_S-  as 
in  equation  5.2.  CS2  and  CS4  are  normalized  by  Sa  as  in  equation  5.4. 


.02-04 

Hz 

.03-05 

Hz 

.04-06 

Hz 

.02-06 

Hz 

CM 

31  (85) 

22  (71) 

22  (64) 

20 (68) 

CS, 

13(41) 

16(47) 

19(53) 

16(45) 

cs2 

17(40) 

19(46) 

22 (52) 

20  (45) 

CS3 

15  (39) 

17(46) 

20  (52) 

17(45) 

cs4 

17(40) 

19(46) 

22  (53) 

20  (45) 

5.3.  The  Effect  of  Signal  Strength  (Event  Size)  on  Backazimuth  Estimates 

Figure  26  shows  the  backazimuth  residuals  for  each  of  the  algorithms  vs.  event  size,  which 
serves  as  a  proxy  for  the  signal-to-noise  ratio.  For  large  events,  the  CM’s  performance  is  almost 
comparable  to  the  CS  algorithm’s,  except  for  a  large  number  of  sign  errors.  The  greatest 
advantage  of  the  CS  algorithm  is  in  estimating  backazimuth  of  smaller  events,  especially  at 
lower  frequency,  where  the  CM  fails  badly. 
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Figure  26.  Median  ±2  SMAD  confidence  intervals  of  backazimuth  residuals  binned  by  Ms 
values  for  the  CM  (dotted)  and  CSi  algorithm  (solid).  Each  bin  has  250  backazimuth  residuals. 
Results  at  0.03-0.05  Hz  are  intermediate  between  those  shown. 

5.4  The  Effect  of  Variable  Group  Velocity  Windows 

The  work  reported  above  used  2.5  to  5.0  km/sec  group  velocity  windows  that  we  expected  to 
favor  the  CM,  as  they  are  most  likely  to  include  both  Love  and  Rayleigh  waves.  Next  we  use  the 
group  velocities  predicted  by  a  1°  global  surface  wave  model  (Stevens,  et  al  2001 )  to  select 
windows  intended  to  more  narrowly  bracket  the  Rayleigh  waves  and  so  minimize  the  effect  of 
noise  outside  those  windows  (Figure  27).  The  minimum  group  velocity  is  that  used  in  the 
dispersion  test  for  surface  waves  at  the  high  frequency  limit  of  the  passband  Similarly,  the 
maximum  group  velocity  is  determined  by  the  low  frequency  limit. 
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Figure  27.  Surface  waves  in  the  3  passbands  tested,  at  the  station  ARCES,  for  an  Ms  5.0  event  at 
6035  km.  The  entire  trace  is  the  5.0  to  2.5  km/s  group  velocity  window.  Outlined  segments  are  the 
shorter  group  velocity  windows  used  to  isolate  the  Rayleigh  waves.  The  isolation  works  best  in  the 
highest  frequency  band,  where  the  large  Love  wave  (top  trace)  is  outside  the  narrower  window. 


Table  4  is  similar  to  Table  3,  but  for  the  variable  length  group  velocity  windows.  The 
performance  of  the  CM  is  slightly  worse  for  the  shorter  windows  at  the  highest  frequency 
passband,  as  expected,  and  there  is  a  small  improvement  in  performance  of  CS|  compared  to  the 
longer  group  velocity  windows.  Overall,  the  effect  of  narrowly  isolating  the  Rayleigh  waves  is 
small,  which  is  encouraging  for  future  implementation  of  the  CS  algorithm  for  routine  detection. 


Table  4.  Same  as  Table  3,  but  for  group  velocity  windows  designed  to  isolate  the  Rayleigh 
waves.  The  change  from  the  2.5  to  5  km/s  group  velocity  windows  is  given  in  parentheses 
(negative  change  is  improvement). 


.02-  04  Hz 

.03-  05  Hz 

.04-  06  Hz 

CM 

31.4  (+0.7) 

21  2  (-0.4) 

23.9  (+1.5) 

CS, 

12.7  (-0.5) 

14.9  (-0.8) 

18.1  (-0.5) 

cs2 

17.0  (+0.5) 

19.0  (+0.3) 

21.5  (-0.3) 

cs3 

14.7  (-0.1) 

16.8  (-0.2) 

19.6  (-0.4) 

cs4 

16.9  (+0.4) 

19.0  (+0.2) 

21.6  (-0.2) 

5.5  Predicting  Accuracy  of  Estimates  Using  the  Cross-Correlation  Value 

The  value  of  the  cross-correlation  of  the  Hilbert  transformed  vertical  with  the  radial,  C?r,  predicts 
the  accuracy  of  the  backazimuth  estimate.  Figure  28  (right)  shows  backazimuth  residuals  vs. 
Cfrfor  CSi,  applied  to  data  that  pass  the  dispersion  test.  The  plot  shows  the  median  ±  two 
SMAD  confidence  intervals  for  bins  with  171  azimuth  residuals  each  versus  the  median 
Cir  value  for  each  bin.  The  histograms  show  azimuth  residuals  for  all  the  data  that  pass  the 

dispersion  test  (left),  and  just  those  with  CIr>  0.8  (right). 
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Figure  28.  Median  azimuth  residuals  ±2  SMAD  confidence  intervals  vs.  the  correlation  of  the 
radial  and  Hilbert  transformed  vertical  Rayleigh  waves  for  the  0.03-0.05  Hz  passband.  Results 
are  similar  for  the  other  passbands. 


At  0.02  to  0.04  Hz,  the  cross-correlations  (CIr)  of  60%  of  the  data  are  >  0  8.  The  SMAD  of  the 
backazimuth  residuals  of  those  data  is  8.5°,  vs.  25.2°  for  data  with  CIr<  0.8.  At  0.03  to  0.05  Hz, 
the  58%  of  the  data  that  have  Cir>  0.8  have  a  SMAD  of  9.3°  vs  28.8°  for  the  data  with  Clr  <  0.8 
At  0.04  to  0.06  Hz,  the  54%  of  the  data  that  have  CIr>  0.8  have  a  SMAD  of  1 1 . 1°  vs  35.8°  for 
the  data  with  Cir<  0.8. 

This  ability  to  predict  backazimuth  estimate  accuracy  can  aid  association  with  known  events.  In 
particular,  the  strict  dispersion  criteria  for  detection  at  the  1DC  could  be  relaxed  in  cases  where  the 
backazimuth  is  consistent  with  the  theoretical  backazimuth  and  the  cross-correlation  value  is  high. 

5.6  Performance  on  Events  for  Which  Suiface  Waves  Are  Not  Detected 

Figure  29  shows  the  results  of  each  method  of  polarization  analysis  applied  to  the  poorer  S/N 
records  that  did  not  pass  the  dispersion  test.  The  CS1  implementation  extracts  accurate 
backazimuth  estimates  for  very  many  of  these  data.  The  more  accurate  estimates  can  be 
identified  by  their  cross-correlation  value. 

The  SMAD  of  the  CS  backazimuth  residuals  for  all  data  not  passing  the  dispersion  test  is  58°,  vs. 
108°  for  the  CM  The  middle  plot  shows  histograms  just  for  data  with  a  cross-correlation  >  0.8, 
which  comprise  14.7%  of  the  data.  Their  backazimuth  residual  has  a  SMAD  of  14.8°.  The  F- 
statistic  of  the  CM  does  not  provide  similar  predictive  capabilities.  In  fact,  for  the  largest  F- 
statistic  values,  the  CM’s  error  increases.  This  is  because  the  F-statistic  can  be  very  large  when 
the  signal  is  dominated  by  very  large  Love  waves,  but  such  records  often  have  1 80°  (i.e.  sign) 
errors.  Even  disregarding  the  sign  errors,  the  F-statistic  is  a  much  poorer  predictor  of  the 
accuracy  of  polarization  estimates,  and  the  resolution  of  the  CM  is  much  poorer  than  that  of  the 
CS  algorithm. 
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Figure  29.  Azimuth  residuals  for  the  0.02  -  0.04  Hz  passband  for  the  CS  algorithm  (top  row) 
and  the  CM  (bottom  row).  Histograms  of  azimuth  residuals  for  all  the  data  (left  column)  and  for 
data  with  correlation  >  0.8  for  CS  (upper  middle)  or  F-statistic  >  30  for  the  CM  (lower  middle). 
The  middle  plots  use  ~14.5%  of  the  data  in  each  case.  The  upper  right  plot  shows  the  median 
azimuth  residual  ±2  SMAD  confidence  intervals  vs.  the  correlation,  as  in  Figure  28,  for  the  CS 
algorithm.  The  bottom  right  shows  the  median  azimuth  residual  ±2  SMAD  confidence  vs.  the  F 
statistic  of  the  CM. 
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6  PERFORMANCE  OF  ALGORITHMS  AND  ANALYSIS  OF 
“PROBLEM  CASES” 

In  many  respects  surface  waves  are  well-behaved  and  straightforward  to  model.  Two 
characteristics  in  particular  are  helpful:  1)  since  the  surface  wave  always  follows  a  minimum 
time  path,  any  path  other  than  the  great  circle  path  is  necessarily  longer,  which  reduces  the  error 
in  the  travel  time  (in  inversions,  the  main  effect  of  this  is  to  make  low  velocity  regions  appear 
somewhat  smaller  than  they  actually  are  (e  g.  Ritzwoller,  et  al,  2002));  and  2)  along  most 
continental  paths  at  the  regional  distances  of  greatest  interest  today,  and  at  the  periods  greater 
than  10  seconds  which  are  recommended  for  discrimination,  amplitude  differences  due  to 
variations  in  attenuation  are  small.  Only  regions  of  unusually  high  attenuation  cause  significant 
errors  in  amplitude  predictions  along  short  paths.  For  most  surface  wave  studies,  uncertainties  in 
event  location  and  measurement  error  are  responsible  for  more  of  the  unmodelable  error  than 
actual  differences  due  to  complexities  in  the  earth.  There  are  some  exceptions  to  this,  however. 

In  cases  where  surface  waves  travel  along  grazing  paths  near  large  changes  in  earth  structure, 
surface  waves  can  be  much  more  complex  and  difficult  to  model.  Multiple  interfering  or  distinct 
surface  wave  arrivals  can  be  observed  that  take  different  paths  through  or  around  a  low  velocity 
region  (Figure  30).  The  Tarim  Basin  in  western  China  is  a  particularly  strong  example  of  this 
effect.  In  this  case,  dispersion  models  can  at  best  predict  the  dominant  surface  wave  arrival,  and 
the  amplitude  may  not  be  predictable  without  detailed  analysis  of  the  data  and  the  earth  structure 
that  causes  the  variations.  One  technique  that  we  found  to  help  for  achieving  consistency  of 
measurement  and  reliability  of  the  inversion  is  to  review  the  data  and  order  (select)  the  arrivals 
such  that  the  dispersion  picks  change  from  fast  to  slow  to  fast  smoothly  as  the  azimuth  changes 
across  the  basin.  Although  this  doesn’t  completely  resolve  the  multiple  arrival  problem,  it  selects 
the  arrivals  in  a  manner  consistent  with  waves  traveling  on  a  continuous  set  of  paths  across  the 
complex  earth  structure.  This  still  does  not  lead  to  a  satisfactory  solution  for  grazing  paths  at 
higher  frequencies,  however,  as  illustrated  in  Figures  31  and  32. 


Figure  30.  Seismograms  recorded  at  two  stations  on  paths  passing  through  the  Tarim  Basin. 
There  are  two  distinct  surface  wave  arrivals  of  comparable  amplitude.  The  second  arrival  is 
dominant  for  paths  traveling  directly  through  the  basin,  and  the  first  arrival  is  dominant  for  paths 
near  the  basin  boundary. 
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Figure  31.  Paths  near  and  across  the  Tarim  Basin.  Colors  correspond  to  the  colors  on  the 
dispersion  curves  in  Figure  32. 


Figure  32.  Observed  (left)  and  calculated  (right)  dispersion  curves  for  the  paths  shown  in  Figure  3 1 . 
Colors  correspond  to  the  paths  on  the  map  in  Figure  3 1 . 


To  determine  the  dispersion  curve  for  a  particular  path,  we  use  the  narrow-band  filtering 
algorithm  originally  due  to  Dziewonski  et  al.  (1969)  as  implemented  by  Stevens  and  McLaughlin 
(2001).  This  method  applies  ~50  narrow  Gaussian  filters  to  a  vertical  seismogram,  calculates  the 
envelope  function,  and  then  determines  the  time  and  instantaneous  frequency  of  each  peak  The 
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peaks  are  due  not  only  to  the  Rayleigh  wave  of  interest,  but  also  to  other  interfering  phases, 
higher  modes,  multi-pathed  arrivals,  and  noise.  We  identify  the  appropriate  dispersion  curve  by 
connecting  the  most  significant  peaks  from  the  set  of  peak  frequencies  and  times,  which  have 
been  converted  to  group  velocities.  It  begins  with  the  largest  peak,  and  then  moves  out  in  both 
directions  finding  the  largest  nearby  peak.  This  works  well  for  cases  where  the  Rayleigh  wave  is 
separated  in  frequency  or  time  from  other  arrivals,  which  is  the  case  for  many  records,  but  will 
not  work  well  over  frequency  bands  in  which  there  are  interfering  arrivals  closely  spaced  in  both 
time  and  frequency.  Situations  such  as  the  Tarim  Basin  can  create  multiple  interfering  surface 
waves  in  some  frequency  bands. 

Figure  33  shows  2  sets  of  3-component  recordings  of  the  same  event,  filtered  from  1 0  to  20 
second  period.  The  propagation  path  (Figure  34)  for  the  record  on  the  left  is  uncomplicated,  and 
the  records  show  a  distinct  Love  wave  on  the  tangential  component,  and  a  distinct  Rayleigh  wave 
on  the  radial  and  vertical  components.  It  appears  that  the  surface  waves  come  in  somewhat  off 
azimuth,  as  there  is  a  little  Love  wave  energy  on  the  radial,  but  there  is  none  clearly  on  the 
vertical  and  the  corresponding  selection  of  a  dispersion  curve  is  straightforward  and  matches  the 
model  well  (Figure  35). 
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Figure  33.  3-component  surface  waves  for  recordings  of  a  single  event  with  one  path  crossing 
(right)  and  one  not  crossing  (left),  the  Tarim  Basin.  Source-station  distances  are  1,427  km  (left) 
and  1,918  km  (right). 
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Figure  34.  Propagation  paths  overlain  on  a  sediment  thickness  map  (contoured  in  kilometers  by 
grayscale  darkness)  for  the  uncomplicated  records  on  the  left  side  of  Figure  33  (propagating  to 
the  northeast)  and  the  complicated  records,  which  propagate  west  through  the  thick  sediments  of 
the  Tarim  Basin. 

The  propagation  path  (Figure  34)  of  the  record  on  the  right  of  Figure  33  traverses  the  Tarim 
Basin  and  appears  to  have  two  Rayleigh  waves  that  interfere  in  the  frequency  band  from  0.04- 
0.06  Hz  and  become  distinct  at  higher  frequencies.  In  this  case,  the  automatic  pick  of  the 
dispersion  curve  jumps  from  higher  to  lower  velocity  between  .045  and  .055  Hz,  as  it  jumps 
from  the  more  coherent  lower  frequency  Rayleigh  wave  to  the  later  (and  larger)  of  the  two  higher 
frequency  Rayleigh  waves  (Figure  35). 
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Figure  35.  Dispersion  curves  (red)  picked  for  the  vertical  component  records  shown  in  Figure 
33,  compared  with  the  model  based  prediction  (maroon). 

In  this  application,  we  are  trying  to  measure  dispersion  curves  to  use  in  tomographic  inversion 
for  earth  structure.  How  do  we  use  data  such  as  that  shown  in  the  right  hand  side  of  Figure  35? 
We  could  simply  throw  out  the  data  as  unreliable,  but  then  would  lose  constraints  on  earth 
structure  in  a  complicated  area.  We  could  decide  to  always  take  the  fastest  arrival,  which 
generally  tends  to  be  more  consistent,  but  that  biases  the  data  set  in  favor  of  arrivals  that  travel 
on  faster  paths,  such  as  around  the  Tarim  Basin  rather  than  through  it.  The  most  disastrous 
approach  is  to  inconsistently  choose  the  faster  or  slower  branches  for  nearby  paths,  in  which  case 
no  consistent  earth  model  is  possible.  The  best  approach  appears  to  be  to  use  the  Rayleigh  wave 
that  travels  closest  to  a  great  circle  path,  and  to  not  use  frequency  bands  where  there  is  strong 
interference.  Following  is  one  method  for  accomplishing  this. 

To  consistently  identify  the  “same”  Rayleigh  wave  from  successive  events,  i.e.  one  that  takes  a 
similar  path,  we  first  identify  the  reason  for  there  being  two  (or  more)  separate  Rayleigh  waves 
in  seismograms  with  paths  crossing  the  Tarim  Basin.  One  hypothesis  is  that  the  later  arrival  has 
taken  the  more  direct,  but  slower  path  through  the  basin,  while  the  earlier  arrival  has  taken  a 
separate  path  around  the  basin.  A  second  possibility  is  that  the  earlier  Rayleigh  wave  is  a 
converted  Love  wave.  That  is,  the  strong  lateral  velocity  contrast  at  the  basin  could  rotate  the 
Love  wave  so  that  it  is  partly  polarized  to  the  vertical.  Similar  observations,  termed  quasi-Love 
waves,  were  made  for  paths  to  Hawaii  (Levin  and  Park,  1998),  in  that  case  attributed  to  lateral 
changes  in  anisotropy. 

The  polarization  algorithm  presented  in  the  previous  section  provides  a  means  of  identifying  the 
Rayleigh  wave  path.  The  correlation  of  the  Hilbert  transformed  radial  with  the  vertical  record 
provides  a  means  of  identifying  Rayleigh  waves  and  estimating  their  polarization  direction. 
Figure  36  (top  row)  shows  the  vertical  seismogram  of  the  problem  case  of  Figure  33,  filtered 
from  10-20  seconds  period  (left)  and  5-10  seconds  period  (right).  In  the  lower  frequency  case, 
the  later  Rayleigh  wave  is  larger.  The  backazimuth  providing  the  maximum  correlation  of  the 
Hilbert  transformed  radial  with  the  vertical  component  is  found  for  a  moving  window  3  times  the 
center  period  of  the  passband.  The  overlap  between  windows  can  be  as  small  as  necessary  for 
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refining  pick  times.  The  maximum  correlation  of  each  window  (second  row)  identifies  both 
Rayleigh  waves,  with  the  second  being  dominant.  The  third  row  shows  the  backazimuth  vs.  time. 
For  the  lower  frequency,  the  backazimuth  of  the  second  Rayleigh  wave,  at  -820  seconds  (2.34 
km/s)  is  88°,  13°  off  azimuth  to  the  north  of  the  theoretical  backazimuth  of  101°.  The 
backazimuth  of  the  earlier  arrival,  at  650  seconds  (2.95  km/s)  is  67°,  34°  off  azimuth  to  the 
north.  This  is  consistent  with  the  earlier  Rayleigh  wave  traveling  around  the  Tarim  Basin  to  the 
north,  while  the  second  one  travels  through  the  basin,  albeit  possibly  also  traveling  through  the 
shallower  part  of  the  basin  somewhat  to  the  north  of  the  great  circle  path. 

The  higher  frequency  records  (Figure  36,  right)  require  a  different  interpretation.  In  that  case, 
the  first  Rayleigh  wave,  at  -650  seconds  (2.95  km/s)  is  much  larger  in  the  seismogram,  and 
corresponds  to  the  large  clear  peak  of  the  correlation  plotted  vs.  time.  That  arrival  however  has  a 
measured  backazimuth  of  165°,  64°  off  azimuth  to  the  south.  Given  the  large  azimuth  deviation 
and  high  velocity,  this  could  be  a  quasi-Love  wave,  that  is  a  Love  wave  refracted  toward  more 
vertical  polarization  by  a  strong  velocity  discontinuity  dipping  perpendicular  to  the  propagation, 
such  as  the  basement  of  the  Tarim  Basin.  Alternately,  the  polarization  could  just  be  due  to 
complexity  of  the  original  Rayleigh  wavefront  due  to  the  basin. 
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Figure  36.  Vertical  component  surface  waves  (top  row)  for  the  problem  path  across  the  Tarim 
Basin.  The  second  row  shows  the  correlation  of  the  Hilbert  transformed  radial  with  the  vertical, 
which  identifies  Rayleigh  waves.  The  third  row  shows  the  backazimuth  vs  time.  The  left  column 
is  for  data  filtered  from  10  to  20  seconds  period,  and  the  right  is  for  data  filtered  from  5  to  10 
seconds  period.  The  great  circle  back  azimuth  is  101  degrees. 
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A  fan  section  of  events  recorded  at  the  KNET  station  ULHL  provides  a  look  at  variations  in  the 
characteristics  of  records  with  paths  through  different  parts  of  the  basin  (Figures  37  and  38). 
Figure  37  shows  the  different  paths  through  the  basin.  We  focus  on  the  events  within  the  circle 
in  the  lower  right. 


Figure  37.  Location  of  paths  through  different  parts  of  the  Tarim  Basin  (red  lines)  from  events 
(within  green  circle  at  lower  right)  to  the  KNET  station  ULHL  (blue  circle).  Sediment  thickness 
in  kilometers  is  indicated  by  grayscale  shading. 


Figure  38  shows  3-component  records,  filtered  from  20  to  30  seconds  period,  of  the  events  with 
paths  through  the  Tarim  Basin,  shown  in  Figure  37.  Records  are  ordered  by  backazimuth,  which 
is  listed  to  the  right  of  the  records.  Traces  are  plotted  on  the  same  scale  for  each  event’s  records, 
normalized  by  the  radial  amplitude.  In  this  passband,  most  of  the  records  look  similar.  The  two 
exceptions  are  the  top  and  bottom  set  of  traces.  The  uppermost  record,  from  the  event  passing 
closest  to  the  center  of  the  basin,  shows  high  amplitude  signal  preceding  the  surface  waves, 
which  may  be  from  an  earlier  event.  The  bottom  set  of  traces  is  for  the  event  coming  closest  to 
the  southern  edge  of  the  basin.  It  has  multiple  peaks,  similar  to  the  problem  case  of  Figure  33, 
which  skirted  the  northern  edge  of  the  basin. 
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20  -  30  Second  Period 


Time  (seconds)  Time  (seconds)  Time  (seconds) 


Figure  38.  3-component  records,  filtered  from  20  to  30  second  period,  arranged  by  backazimuth 
of  the  events  within  the  green  circle  in  the  lower  right  of  Figure  37.  Records  are  centered  at  2.8 
km/sec  (0  seconds  on  the  time  axis).  Fifty  seconds  either  side  of  the  center  are  highlighted  in 
each  record  to  facilitate  comparison. 

10-20  Second  Period 


Time  (seconds)  Time  (seconds)  Time  (seconds) 


Figure  39.  Same  as  Figure  38,  except  filtered  from  10  to  20  second  period. 
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Figure  39  shows  the  same  records  filtered  from  10  to  20  second  period.  The  influence  of  the 
basin  appears  to  be  stronger,  as  the  records  vary  much  more.  In  some,  for  example  at  125° 
backazimuth,  it  appears  that  there  could  be  quasi-Love  waves  on  the  vertical  component.  Many 
of  the  records  appear  to  be  complicated  by  multipathing.  The  differences  between  records 
however,  are  not  simple,  systematic,  or  easily  predictable.  In  the  following  we  use  Rayleigh 
wave  polarization  to  try  to  resolve  the  different  arrivals. 

Figure  40  (left  column)  shows  the  vertical  and  Hilbert  transformed  radial  records,  filtered  from 
20  to  30  seconds  period,  of  Figures  37-39.  The  right  column  shows  their  crosscorrelation.  This 
differs  from  the  example  in  Figure  36,  which  showed  the  maximum  crosscorrelation  in  each 
time  window  found  at  different  backazimuths  for  each  window.  Since  the  crosscorrelation  is 
formed  using  the  radial  component  rotated  to  the  great  circle  path,  this  example  is  biased  toward 
identifying  Rayleigh  waves  arriving  at  the  predicted  backazimuth.  In  this  frequency  band,  there 
is  little  ambiguity  about  the  phase  selection  and  arrival  time  for  most  paths.  However,  in  the 
126.9°  case  (bottom  row),  there  is  an  early,  possible  quasi-Love  wave  arrival  and  a  later, 
possibly  multi-pathed  arrival  in  the  vertical  seismogram,  but  there  is  only  one  distinct  arrival  in 
the  correlation. 
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Figure  40.  Vertical  (blue)  and  Hilbert  transformed  radial  (dotted  red)  seismograms  (left) 
filtered  from  20  to  30  seconds  period,  and  their  correlation  over  a  series  of  tightly  overlapping  75 
second  windows. 
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Figure  41  is  similar  to  Figure  40,  but  for  10  to  20  seconds  period.  These  records  are  more 
complex  than  those  at  20  to  30  seconds  period  in  Figure  6.11,  and  the  correlations  of  the  Hilbert 
transformed  radial  and  the  vertical  seismograms  are  also  more  complex.  Starting  this  time  at  the 
southern  end,  the  vertical  component  record  at  126.9°  (bottom  row),  has  3  distinct  arrivals,  but 
the  correlation  shows  only  two  are  consistent  with  a  Rayleigh  wave  arriving  at  or  near  the 
predicted  backazimuth.  The  next  4  records  to  the  north  also  have  two  distinct  arrivals  in  the 
correlation  vs.  time.  Those  at  124  1°  and  124.4°  backazimuth  have  some  additional  complexity, 
but  less  than  the  seismograms.  The  seismograms  also  appear  to  have  very  early  quasi-Love 
wave  arrivals,  which  are  diminished  in  the  crosscorrelations.  To  distinguish  between  the  two 
Rayleigh  waves,  and  associated  peaks  from  one  record  to  the  next,  we  need  to  do  more.  We 
could  apply  the  algorithm  used  to  make  in  Figure  36,  to  estimate  the  backazimuth  associated 
with  each  Rayleigh  wave  coming  in  near  the  predicted  backazimuth  and  directly  estimate  the 
dispersion  curve  from  the  correlations  Alternately,  we  could  estimate  the  backazimuth  of  each 
arrival  peak  shown  in  Figure  35.  We  could  then  eliminate  from  consideration  any  arrivals  not 
arriving  from  a  direction  near  the  predicted  backazimuth,  and  use  just  the  remaining  peaks. 


Vertical  and  HRadial  10-20  Second  Period  Crosscorrelation 
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Figure  41.  Vertical  (blue)  and  Hilbert  transformed  radial  (dotted  red)  seismograms  (left) 
filtered  from  10  to  20  seconds  period,  and  their  correlation  over  a  series  of  tightly  overlapping  45 
second  windows. 
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We  can  examine  the  utility  of  such  approaches  by  determining  the  backazimuth  that  provides  the 
maximum  crosscorrelation  of  the  vertical  and  Hilbert  transformed  radial  vs.  time,  as  in  Figure  36, 
for  the  peaks  of  the  two  southernmost  records  of  Figure  41,  which  are  at  1,771  km  (126.9° 
backazimuth  event)  and  1457  km  (125.0°  backazimuth  event)  respectively.  Figure  42  shows 
those  results.  These  crosscorrelation  plots  represent  the  crosscorrelation  at  the  rotation  of  the 
horizontal  components  that  maximizes  the  crosscorrelation.  Those  of  Figure  41  are  the 
crosscorrelation  using  just  the  radial  component  given  the  theoretical  backazimuth.  The  first 
peak  of  the  southernmost  records’  crosscorrelation  vs.  time  from  the  origin  arrives  at  2.78  km/s 
(637  seconds)  and  is  47°  off  azimuth  to  the  south.  Thus  it  is  either  refracted  strongly  off  the  great 
circle  path,  or  is  a  quasi-Love  wave,  or  perhaps  both  such  arrivals  overlap  The  second,  at  2.48 
km/s  (715  seconds)  is  just  4°  off  azimuth  to  the  south,  and  so  represents  the  direct  Rayleigh 
wave.  The  same  analysis  on  the  next  event  to  the  north,  at  125°  backazimuth,  finds  that  the  first 
peak  of  the  crosscorrelation,  at  2.66  km/s  (548  seconds)  is  45°  off  azimuth  to  the  south.  Thus  it 
appears  to  correspond  with  the  first  arrival  of  the  next  event  to  the  south.  The  second  arrives  at 
2.26  km/s  (645  seconds)  and  is  23°  off  azimuth  to  the  south.  This  event’s  path  is  through  a 
deeper  part  of  the  basin,  so  it  has  slower  propagation  and  greater  deflection  (toward  the 
shallower  sediments  further  south).  This  demonstrates  that  polarization  can  be  used  to  identify 
and  associate  distinct  arrivals  between  events,  and  should  improve  the  determination  of  accurate 
and  appropriate  group  velocity  curves. 
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Figure  42.  Vertical  component  surface  waves  (top  row)  for  the  two  southernmost  records  of 
Figure  41  (126.9°  backazimuth  on  the  left,  125.0°  backazimuth  on  the  right),  the  correlation  of 
the  horizontal  component  rotated  to  the  backazimuth  that  maximizes  the  crosscorrelation  with 
the  Hilbert  transformed  vertical  for  successive  windows  (middle  row),  and  that  backazimuth  vs 
time  (bottom  row).  The  horizontal  line  indicates  the  theoretical  backazimuth.  Data  are  filtered 
from  10  to  20  seconds  period.  Time  is  seconds  from  origin. 


50 


7  CONCLUSIONS  AND  RECOMMENDATIONS 


The  dispersion  models  developed  in  this  study  can  be  used  to  find  the  phase  and  group  velocities 
between  any  two  points  on  the  earth  at  any  frequency.  These  models  are  available  to  all 
researchers  on  request,  and  have  been  provided  for  use  in  a  number  of  other  research  projects. 
The  dispersion  models  depend  on  dispersion  data  for  accuracy,  and  so  can  be  improved  as  more 
dispersion  data  becomes  available.  We  therefore  recommend  continuation  of  this  iterative 
process.  As  more  data  becomes  available,  particularly  from  regions  or  in  frequency  bands 
currently  sparsely  covered,  we  can  improve  the  models  and  in  turn  provide  these  predictions  for 
comparison  with  and  as  a  starting  point  for  new  surface  wave  studies  in  regions  of  interest. 

The  azimuth  estimation  algorithm  currently  used  at  the  International  Data  Center  should  be 
replaced  with  the  Selby/Chael  algorithm  that  was  implemented  and  tested  in  this  project.  The 
procedure  can  be  used  together  with  the  dispersion  test  to  reduce  the  detection  threshold  for 
automatic  surface  wave  processing  by  relaxing  the  requirements  of  the  dispersion  test  when  the 
azimuth  estimate  is  consistent  with  an  event  hypothesis  with  a  high  correlation  value. 

The  path  corrected  spectral  magnitude  provides  a  “regional  Ms”  which  has  been  a  long-term  goal 
of  this  program.  The  path  corrected  spectral  magnitude  can  be  measured  over  any  frequency 
band,  however  for  discrimination  purposes,  surface  waves  should  be  measured  at  periods  greater 
than  10  seconds  Periods  of  10  seconds  and  longer  can  be  measured  even  at  very  close  range,  and 
there  is  only  a  small  S/N  improvement,  if  any,  at  shorter  periods.  Measurement  at  higher 
frequencies  risks  error  in  discrimination  due  to  the  decrease  in  earthquake  spectral  amplitude 
relative  to  explosions  at  higher  frequencies,  particularly  for  events  for  which  shallow  depth  is 
uncertain.  A  path  corrected  time  domain  magnitude  can  be  derived  by  combining  the  path 
corrected  spectral  magnitude  with  the  Russell  (2004)  surface  wave  magnitude  M«(b). 

Additional  work  is  needed  regarding  how  to  model  and  invert  surface  wave  amplitude  and 
dispersion  in  regions  of  complex  earth  structure.  In  particular,  surface  waves  traveling  at  near 
grazing  incidence  to  strong  heterogeneities  may  exhibit  amplitude  variations  not  only  due  to 
focusing/defocusing,  but  also  due  to  interference  of  arrivals  taking  multiple  paths  through  the 
medium.  Similarly,  dispersion  measurements  are  complicated  by  the  presence  of  distinct, 
multiple  arrivals  at  short  periods  and  interfering  arrivals  that  cannot  be  separated  at  longer 
periods.  Some  of  this  behavior  can  be  modeled  by  incorporating  scattering  theory  into  the 
inversion  procedure  (e  g.  Zhou  et  al.,  2004),  however  the  interference  is  so  strong  in  some  cases 
that  it  may  not  be  possible  to  model  it  accurately  with  approximate  methods.  A  partial  solution  is 
to  use  waveform  polarization  to  identify  waves  traveling  closest  to  a  great  circle  path  and  to  use 
only  frequency  bands  where  polarization  is  reasonably  consistent  with  a  great  circle  path  and 
interference  between  multiple  arrivals  is  not  present  or  can  be  removed. 
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DATA  DELIVERABLE 


Two  data  deliverables  are  provided  together  with  this  final  report.  They  are: 

1 .  The  final  set  of  earth  models,  dispersion  curves  and  other  derived  data.  These  are 
contained  in  compressed  tar  file  “LP_2005_Mar.tar.Z’\  Information  on  data  format  is 
included  in  the  delivery. 

2.  The  maxpmf  program  compiled  for  Sun  Solaris,  and  the  maxpmf  man  page.  These  are 
contained  in  compressed  tar  file  “maxpmf_5.3.tar.Z” 

These  data  deliverables  may  be  obtained  from  the  contracting  office  or  from  the  authors. 
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